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ABSTRACT 


A  two-level,  quasi-geos trophic  atmospheric  model  on  a 
mid-latitude  $-plane  with  a  single  wave  in  the  x-direction 
was  coupled  with  a  three  level  primitive  equation  ocean 
model  through  the  exchange  of  sensible  heat  and  momentum 
to  investigate  the  effect  of  the  ocean  circulation  on  tran- 
sients in  the  atmospheric  circulation  and  vice  versa.   A 
long-term  (79  years)  simulation  of  the  atmosphere-ocean 
system  was  carried  out.   The  experiment  was  conducted  in 
three  phases  that  consisted  of  integrating  the  combined, 
the  ocean-only  and  the  atmosphere-only  model  equations. 

The  combined  model  produced  realistic  patterns  of  both 
atmospheric  and  oceanic  circulations.   Energy  studies  of 
the  atmosphere  were  completed  with  the  results  comparing 
favorably  with  the  existing  knowledge  of  the  energy  balance 
in  the  real  atmosphere. 

Data  were  retained  at  selected  points  on  a  daily  basis 
and  spectrally  analyzed  for  a  simulated  six  year  period. 
Transients  with  periods  on  the  order  of  one  and  two  months 
appeared  in  the  spectra  of  both  atmospheric  and  oceanic 
variables  from  the  combined  model.   These  transients  were 
not  observed  in  the  ocean-only  or  atmosphere-only  phases  of 
the  experiment,  and  are  attributed  to  the  interaction  occur 
ring  in  the  combined  model. 
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I.   INTRODUCTION 


The  use  of  numerical  models  to  simulate  physical  con- 
cepts that  would  otherwise  be  restricted  to  limited  observa- 
tional studies  has  advanced  the  knowledge  of  both  the 
atmosphere  and  ocean  and  their  large  scale  interaction  a 
great  deal  over  the  last  several  years.   But  that  knowledge 
is  far  from  complete  and  it  is  the  intent  of  this  work  to 
reveal  more  information  on  the  long  term,  large  scale  inter- 
actions and  the  frequencies  of  the  response  between  the 
atmosphere  and  ocean.   The  primary  goal  of  this  study  was 
to  identify  the  effects  of  the  ocean  circulation  on  the 
transient  (time-dependent)  atmospheric  circulation,  and 
similarly  the  effect  of  the  atmospheric  circulation  on  the 
ocean . 

The  development  of  a  large  scale  air-sea  interaction 
model  that  would  remain  computationally  stable  throughout 
a  long-term  numerical  integration  of  the  dynamical  equa- 
tions of  motion  was  required  to  achieve  this  goal.   It  was 
desired  that  the  method  by  which  the  integrations  proceeded 
be  one  which  allowed  the  synchronization  of  the  oceanic  and 
atmospheric  time  steps  so  that  a  "one-to-one"  relationship 
existed  between  the  two  media.   Many  complex  numerical 
models  of  both  the  ocean  and  the  atmosphere  have  been 
developed  in  the  past.   However,  those  models  all  require 
large  amounts  of  computer  time  because  of  the  very  short 


time  steps  required  to  satisfy  the  linear  computational  sta- 
bility requirement.   The  model  developed  in  this  work  had 
fewer  degrees  of  freedom  in  the  atmosphere  and  this  resulted 
in  faster  integration  of  the  atmospheric  equations.   It  also 
utilized  certain  approximations  and  assumptions  that  permit- 
ted the  use  of  longer  time  increments  in  the  integrations. 
The  necessity  of  employing  a  technique  to  simulate  time 
synchronization  was  thereby  avoided. 

Many  numerical  general  circulation  models  have  been 
developed  in  the  last  two  decades.   Phillips  (1956)  made  a 
long  period  forecast  with  a  two-level,  quasi-geos trophic 
model,  starting  with  an  atmosphere  at  rest.   Diabatic  ef- 
fects were  included  in  the  equations  as  a  linear  function 
of  latitude,  and  such  quantities  as  the  intensity  of  the 
heating,  the  value  of  the  vertical  stability  and  the  type  of 
frictional  dissipation,  were  empirically  specified.   The 
flow  patterns  which  developed  were  quite  realistic.   With 
the  development  of  this  model,  Phillips  demonstrated  the 
feasibility  of  numerical  simulation  of  the  atmospheric 
general  circulation  with  a  simple  model. 

Mintz  (1958)  proposed  numerical  general  circulation 
experiments  that  specified  a  heating  function  in  which  the 
diabatic  heating  of  the  air  depended  on  the  air-sea  tempera- 
ture difference.   He  showed  that  this  heating  function  could 
be  used  to  study  the  circulation  of  the  atmosphere  on  a 
planet  covered  by  either  land  or  ocean  or  a  combination 
thereof . 


In  both  of  the  above  works,  the  partial  differential 
equations  were  formulated  in  finite  differences.   Lorenz 
(1960)  discussed  a  method  of  expressing  a  dependent  variable 
as  a  series  of  orthogonal  functions.   An  example  of  such  an 
analysis,  in  one  dimension,  is  the  Fourier  series,  in  which 
the  orthogonal  functions  are  sines  and  cosines.   The  coeffi- 
cients of  these  orthogonal  functions  become  the  new  dependent 
variables  and  a  finite  number  of  variables  are  retained  for 
prediction.   Bryan  (1959)  reported  on  his  work  with  Lorenz 
in  which  this  technique  was  applied  to  the  variables  in  a 
two-layer,  quasi-geos trophic  model.   The  variables  were  ex- 
panded into  a  few  leading  terms  of  a  spherical  harmonic 
series  to  conduct  a  simple  general  circulation  experiment. 

Smagorinsky  (1963)  reported  on  an  extended  period  numeri- 
cal integration  of  a  two-level,  baroclinic  primitive  equa- 
tion model  which  had  "land"  as  its  underlying  surface.   This 
experiment  was  made  for  the  simulation  and  study  of  the 
dynamics  of  the  atmospheric  general  circulation.   The  solu- 
tions corresponding  to  external  gravitational  oscillations 
were  filtered  by  requiring  the  vertically  integrated  diver- 
gence  to  vanish  identically.   This  was  a  major  departure 
from  earlier  work  in  that  many  hydrodynamic  and  kinematic 
constraints  were  removed. 

Mintz  (1964)  conducted  an  experiment  in  climate  simula- 
tion by  numerically  integrating  the  primitive  equations  of 
atmospheric  motion,  with  friction  and  heating  terms  included, 
over  a  long  period  of  time.   His  two-level  model  domain 


encompassed  the  entire  earth,  with  the  underlying  oceans 
and  continents  having  different  thermal  properties.   Leith 
(1965)  described  a  six-level  numerical  model  with  which  an 
attempt  was  made  to  compute  the  hydrodynamic  and  thermo- 
dynamic evolution  of  a  moist  atmosphere  on  the  whole  globe. 
The  model  took  into  account  such  external  influences  as 
solar  heating,  evaporation  and  surface  friction. 

Smagorinsky,  Manabe  and  Holloway  (1965)  reported  the 
results  of  a  nine-level  general  circulation  model  of  the 
atmosphere.   The  primitive  equations  were  used  in  that 
study  in  which  the  atmosphere  attained  a  state  of  quasi- 
equilibrium.   The  nine  levels  of  the  model  were  distributed 
to  resolve  surface  boundary  layer  fluxes  as  well  as  radia- 
tive transfer  by  the  prescribed  distribution  of  ozone, 
carbon  dioxide  and  water  vapor.   This  model  was  an  outgrowth 
of  Phillips'  two-level  quasi-geos trophic  model  and  Smagor- 
insky's  two-level  primitive  equation  model. 

Arakawa  (1966)  derived  a  finite  difference  Jacobian 
expression  for  the  advection  term  in  the  difference  vorti- 
city  equation  for  two-dimensional  incompressible  flow  that 
properly  represented  the  interaction  between  grid  points. 
He  showed  that  the  use  of  this  differencing  technique  main- 
tained the  integral  constraints  on  quadratic  quantities  of 
physical  importance  and  that  this  was  sufficient  to  prevent 
nonlinear  computational  instability.   Before  Arakawa's  work, 
long-term  integrations  had  been  possible  only  with  the  use 
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of  large  friction  and  diffusion  terms.   Arakawa ' s  scheme 
permitted  long-term  integrations  without  the  requirement  for 
such  artificialities. 

Kasahara  and  Washington  (1967)  described  a  model  of  the 
general  circulation  of  the  earth's  atmosphere  which  was 
developed  at  the  National  Center  for  Atmospheric  Research 
(NCAR) .   A  distinguishing  feature  of  this  model  was  that 
the  vertical  coordinate  was  height  rather  than  pressure, 
though  hydrostatic  equilibrium  was  maintained  in  the  system. 
Experiments  with  a  two-level  version  of  this  model  were  con- 
ducted.  Long-term  (over  100  days)  integrations  were  conduc- 
ted testing  various  physical  parameters  and  examining  their 
affect  on  the  results. 

Numerical  models  have  also  been  developed  for  predicting 
the  atmospheric  conditions  in  enough  detail  to  provide 
daily  forecasts  for  given  areas.   Kesel  and  Winninghoff 
(1972)  reported  on  the  U .  S.  Navy  operational  five-level 
primitive  equation  model  that  has  been  used  successfully 
for  operational  forecasts  since  1970.   This  hemispheric 
model  perhaps  typifies  the  state  of  sophistication  to  which 
numerical  weather  prediction  has  progressed. 

Early  mathematical  models  of  the  ocean  circulation  were 
very  simple  as  were  the  early  atmospheric  models.   Munk 
(1950)  utilized  the  oceanic  mass  transport  to  examine  a 
wind-driven  ocean  which  did  not  include  specification  of 
the  vertical  distribution  of  density  or  currents.   He  ad- 
hered to  the  linear  theory  in  obtaining  analytic  solutions 
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to  the  equations  for  steady  horizontal  motions.  Fofonoff 
(1954)  studied  f ric tionless ,  steady  horizontal  flow  which 
could  occur  in  a  homogeneous  ocean  of  constant  depth.  The 
effects  of  the  thermohaline  structure  and  compressibility 
were  neglected  everywhere.  By  considering  this  idealized 
case  he  avoided  the  analytical  difficulties  presented  by 
the  nonlinear  terms  in  the  hydr odynami cal  equations. 

Veronis  (1963a,  1963b,  1965,  1966)  conducted  numerous 
numerical  experiments  using  both  a  grid  point  method  of 
examination  and  a  Fourier  method  of  analysis  and  made  com- 
parisons of  the  two  methods.   He  generated  flow  patterns  in 
a  rectangular,  3-plane,  barotropic  ocean  by  assuming  the 
basin  could  be  divided  into  an  interior  region  in  which  the 
dynamical  balance  was  quas  i*geos trophic  and  a  boundary  region 
where  the  dynamical  balance  was  governed  by  inertial  forces 
much  the  same  as  Munk .   He  developed  both  linear  and  non- 
linear models  and  experimentally  determined  optimum  values 
of  parameters,  such  as  the  Rossby  number  and  the  bottom 
friction  parameter  needed  for  the  numerical  work. 

Bryan  (1963)  obtained  solutions  for  a  time  dependent, 
nonlinear,  barotropic  model  of  a  wind  driven  ocean  by 
numerical  integration.   He  varied  parameters  such  as  the 
wind  stress,  the  Rossby  number  for  the  interior  flow  and 
an  effective  Reynolds  number  for  the  boundary  current  and 
examined  the  results.   He  examined  boundary  currents  and 
the  separation  of  the  western  boundary  current  from  the 
coast  for  both  linear  and  nonlinear  cases.   Bryan  and  Cox 
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(1967)  investigated  the  ocean  circulation  using  a  baro- 
clinic  model  to  numerically  integrate  the  primitive  equa- 
tions in  almost  complete  form  for  different  values  of 
Reynolds  numbers.   The  near  linear  solutions  for  low  Rey- 
nolds number  were  most  accessible  to  a  numerical  net  with 
coarse  resolution  and  gave  valuable  information  for  under- 
standing the  more  complicated  regimes  at  higher  Reynolds 
number . 

Gates  (1968)  integrated  the  primitive  equations  for  a 
rectangular  homogeneous  ocean  with  a  free  surface  on  a  $- 
plane  starting  from  an  initial  state  of  rest.   He  investiga- 
ted the  transient  Rossby  waves  and  found  a  well-marked  life 
cycle  of  amplification  and  decay  as  they  propagated  west- 
ward  at  about  1  m  sec 

Bryan  (1969)  continued  his  earlier  work  with  Cox  by 
adding  new  features  to  their  five-level  model.   He  treated 
salinity  separately  from  temperature  and  allowed  sea  ice 
formation  in  the  new  model.   This  work  is  discussed  in  more 
detail  later  in  the  section  on  the  combined  model.   Haney 
(1971)  discussed  the  formulation  of  the  net  downward  heat 
flux  at  the  ocean  surface,  and  used  this  formulation  in 
experiments  with  a  six-level  baroclinic  ocean  model  similar 
to  Bryan's  with  emphasis  on  the  equatorial  regions.   This 
model  is  representative  of  the  present  state  of  knowledge 
in  ocean  modeling. 

In  recent  years  research  on  the  problem  of  atmosphere- 
ocean  interaction  has  progressed  to  the  point  where  very 
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complex  atmospheric  numerical  models  have  been  combined  with 
similarly  complex  ocean  models  and  integrated  for  extended 
periods  of  time.   The  major  work  in  this  area  combined  a 
multilevel,  primitive  equation  atmospheric  model  including 
the  earth's  hydrology  reported  on  by  Manabe  (1969)  with  a 
multilevel,  primitive  equation  ocean  model  (Bryan  1969). 
The  goal  of  their  study,  quite  different  from  the  present 
one,  was  to  determine  the  effect  of  the  ocean  circulation 
on  the  time  average  climate  of  the  atmosphere.   Their  ex- 
periments produced  some  notable  results  such  as  realistic 
heat  transports  and  reasonable  distribution  of  the  pole- 
ward transport  of  energy.   In  one  sense,  however,  this  ex- 
periment was  considered  incomplete  in  that  the  goal  of 
calculating  a  climate  equilibrium  with  the  combined  system 
was  not  completely  attained. 

That  research  was  carried  out  in  three  steps: 

a.  A  study  of  the  atmospheric  model  without  the  effect 
of  ocean  circulation.   The  effect  of  heat  transfer  by  the 
ocean  current  was  neglected.   The  ocean  had  zero  heat 
capacity  and  simply  acted  as  an  infinite  reservoir  of 
moisture  for  the  atmosphere. 

b.  A  study  of  the  oceanic  model  without  any  feedback 
from  the  atmosphere.   The  distribution  of  surface  tempera- 
ture, wind  stress  and  precipitation  over  the  ocean  were 
given  and  held  constant  with  respect  to  time. 
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c.   A  study  of  the  joint  ocean-atmosphere  model  in  which 
the  two  systems  were  allowed  to  interact  fully  with  each 
other . 

A  comparison  of  the  results  of  the  three  experiments 
yielded  some  insight  into  the  role  of  the  ocean  in  maintain- 
ing the  climate  and  of  the  atmosphere  in  maintaining  the 
thermal  and  dynamical  structure  of  the  ocean. 

In  the  Manabe-Bryan  model  more  computing  time  was  needed 
to  obtain  a  one  day  integration  with  the  atmospheric  model 
than  with  the  ocean  model.   Therefore,  attempts  were  made 
to  economize  computer  time  by  integrating  the  atmospheric 
equations  for  only  one  year  while  integrating  the  oceanic 
equations  for  100  years.   This  resulted  in  a  distorted 
response  of  the  ocean  circulation  to  fluctuations  of  atmos- 
pheric circulation,  since  an  atmospheric  disturbance  that 
should  have  an  effect  for  possibly  one  week  influenced  ocean 
conditions  for  up  to  two  years.   For  example,  it  took  one 
week  for  a  cyclone  wave  to  traverse  the  oceanic  region. 
Because  of  the  synchronization,  the  oceanic  part  of  the 
model  felt  that  cyclone  for  as  long  as  two  oceanic  model 
years . 

Wetherald  and  Manabe  (1972)  reported  on  the  effect  of 
seasonally  varying  the  solar  radiation  in  the  combined 
model  just  discussed.   Many  of  the  results  obtained  in  the 
new  version  of  the  model  were  different,  the  most  signifi- 
cant of  which  was  a  warming  trend  which  occurred  when 
seasonal  variations  of  the  solar  radiation  were  introduced. 
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It  should  be  pointed  out  that  an  exact  statistical  equili- 
brium was  not  obtained  in  any  part  of  the  region  of  this 
integration . 

After  an  extensive  review  of  the  literature  which  in- 
cluded all  of  the  above  works,  it  was  decided  that  to  reach 
a  state  of  quasi-equilibr ium  with  a  coupled  and  time- 
synchronous  atmosphere-ocean  model  the  simplest  atmosphere- 
ocean  system  which  still  contained  the  essential  features 
of  the  dynamics  of  each  media  was  required.   The  atmospheric 
model  selected  was  a  version  of  Phillips'  (1956)  two  level, 
quas i-geos trophic  model  modified  by  implementing  Lorenz's 
(1960)  idea  of  replacing  the  dependent  variable  with  Fourier 
coefficients.   This  decision  was  based  primarily  on  the 
speed  with  which  the  equations  of  this  model  could  be 
integrated.   The  equations  governing  the  atmosphere  could 
be  integrated  for  one  model  year  in  about  30  minutes  on  the 
IBM  360/67.   A  model  similar  to  the  one  desired  had  been 
used  in  an  earlier  atmosphere-only  study  at  the  Naval  Post- 
graduate School  by  Taylor  (1970)  and  was  readily  available 
for  further  modification.   The  details  of  the  atmospheric 
model  are  described  in  Section  II. 

The  ocean  model  selected  was  a  modified  form  of  Haney's 
(1971a)  six  level  baroclinic  model.   The  number  of  levels 
was  reduced  to  three  and  the  horizontal  domain  reduced  to 
an  area  in  the  Northern  Hemisphere.   The  depth  was  reduced 
and  it  was  further  simplified  by  using  the  $-plane  approxi- 
mation.  While  this  primitive  equation  model  was  fairly 
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complex,  it  was  readily  available  for  modification.   De- 
tails of  the  ocean  model  are  described  in  Section  III. 

The  ocean  and  atmosphere  models  were  joined  through 
the  exchange  of  momentum  and  sensible  heat.   Unlike  the 
model  of  Manabe  (1969a),  the  atmosphere  contained  no  mois- 
ture so  there  was  no  explicit  exchange  of  moisture.   In- 
stead, the  effects  of  evaporation  and  precipitation 
processes  were  simulated  by  other  means.   These  simulations 
as  well  as  the  details  of  the  coupling  of  the  two  models 
are  described  in  Section  IV. 

In  a  manner  somewhat  similar  to  the  experiment  of 
Manabe  and  Bryan,  the  research  was  conducted  in  three  phases. 
In  phase  I,  equations  of  the  combined  model  were  integrated 
until  a  state  of  quasi-equilibrium  was  achieved.   This 
required  79  years  of  simulation  using  more  than  100  hours  on 
the  IBM  360/67  computer.   The  solution  of  a  complex  non- 
linear system  of  equations  by  numerical  integration  over  a 
long  time  period  as  was  conducted  in  this  simulation  caused 
some  concern  over  the  accuracy  of  the  numerical  results. 
Error  analysis  studies  by  Fisher  (1965)  and  others  have 
shown  that,  even  for  short  time  integrations,  errors  in 
amplitude  and  phase  arise  in  the  numerical  solution.   This 
is  not  to  say  that  such  errors  make  the  results  meaningless 
because  even  if  fairly  large  errors  were  present  at  the 
end  of  a  long-term  simulation  of  the  atmospheric  circulation, 
the  gross  features  of  the  solution  would  not  be  changed 
significantly.   However,  the  results  should  be  examined  more 
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from  a  qualitative  than  a  quantitative  point  of  view.   The 
results  in  this  long-terra  experiment  surely  contain  numeri- 
cal errors  but  these  results  are  considered  to  be  satisfac- 
tory for  this  study  of  the  coupled  atmosphere-ocean  system. 

One  hundred  day  averages  of  the  atmospheric  surface 
stress  and  surface  temperature  were  calculated  and  used  as 
driving  mechanisms  for  the  ocean  in  phase  II,  an  ocean-only 
experiment.   The  time  of  this  experiment  corresponded  to  the 
last  six  years  of  phase  I.   Similarly  averaged  sea  surface 
temperatures  were  used  as  a  forcing  function  during  phase 
III,  an  atmosphere-only  experiment.   This  also  was  initial- 
ized so  that  the  time  corresponded  to  the  last  six  years  of 
phase  I.   Details  of  the  procedures  used  in  all  three 
phases  are  described  in  Section  V. 

The  resulting  data  were  examined  in  two  ways.   The  first 
was  a  synoptic- type  study  in  which  maps  of  both  ocean  and 
atmospheric  parameters  were  studied  and  their  interactions 
noted.   The  combined  model  produced  both  atmospheric  and 
oceanic  features  that  compared  favorably  with  observed  fea- 
tures.  Energy  balance  diagrams  were  constructed  for  the 
atmosphere  and  these  compared  well  with  existing  knowledge 
of  the  real  atmosphere.  These  synoptic  results  are  shown 
in  Section  VI  A. 

The  second  method  of  examining  results  was  a  spectral 
method.   Time  series  of  data  recorded  over  the  last  six 
years  for  selected  locations  in  the  ocean  and  atmosphere 
were  spectrally  analyzed  to  determine  if  any  short  or 
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intermediate  term  fluctuations  such  as  those  studied  by 
Bjerknes  (1962)  and  Namias  (1969)  appear  due  to  the  inter- 
actions.  The  results  of  this  analysis  are  discussed  in 
Section  VI  B  of  this  work.   By  comparing  the  spectral 
characteristics  of  the  coupled  atmosphere  with  those  of 
the  atmosphere  with  fixed  ocean  surface   temperatures,  the 
role  of  the  ocean  circulation  in  producing  certain  tran- 
sient atmospheric  responses  may  be  studied. 
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II.   THE  ATMOSPHERIC  MODEL 

A.   DESCRIPTION 

The  atmospheric  model  was  based  on  the  selective  appli- 
cation of  the  quasi-geos trophic  theory  to  the  equations  of 
motion.   Since  the  large-scale  motions  in  extra- tropical 
latitudes  are  quasi-geos trophic  in  nature,  the  domain  of 
this  model  was  limited  to  the  middle  latitudes.   The  model 
was  constructed  on  a  3-plane  with  a  central  latitude  of  45 
degrees.   Phillips  (1963)  discussed  geostrophic  motions 
that  were  characterized  by  approximate  uniformity  of  poten- 
tial vorticity.   These  motions  have  characteristic  wave- 
lengths of  about  4000  km.   He  used  an  expansion  in  the 
Rossby  number  about  latitude  45  degrees  to  show  that  the 
3-plane  approximation  is  justified  for  geostrophic  motions 
of  this  type . 

Phillips  (1956)  developed  a  two-level  quas i-geos trophic 
model  that  he  considered  was  the  simplest  model  in  which 
heating  and  friction  could  be  incorporated.   The  basic  model 
described  in  this  thesis  followed  Phillips'  work  in  the 
manner  in  which  the  quasi-geos trophic  equations  were  derived. 

The  quasi-geos trophic  vorticity  equation  on  a  (3-plane  ±s 

1^    (V2i|>)    +  IkxVij;    •    V(V2^)    +   P   ||  -    fQ  ~ 

=  -  8  4-  It-2-  -  a^J.  (2-1) 

dp       dx  dy 
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where  T   is  the  component  of  the  stress  in  the  j\   direction, 
x 

T   is  the  component  of  the  stress  in  the  jj  direction 
and  other  symbols  have  their  usual  meteorological  meaning. 

The  model  is  restricted  to  two  levels  with  a  vertical 
structure  as  shown  in  Figure  1.   The  stream  function  ^ ..  and 
i|;_  are  representative  of  the  flow  patterns  in  the  upper  and 
lower  layers  of  the  atmosphere  respectively. 


■250 


•500 


-750 


■1000 


Figure  1.   Vertical  structure  of  atmospheric  model. 


and 


Boundary  conditions  at  the  top  and  bottom  are: 

a)  =  0   at   p  =  0 
03  =  0   at   p  =  p ,  . 


(2.2) 


The  stress  IT   vanishes  at  p  =  0 .   The  stress  at  anemometer 


level  was  computed  from 

,T  =  P(CD)4  |\V|  W 


(2.3) 
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where  (C  ),  is  a  constant  drag  coefficient  and  \V  is  the  wind 
at  anemometer  level.   The  stress  at  level  four  was  obtained 
by  assuming  that  \V  was  a  fraction,  J(  >  of  the  vectorially  ex- 
trapolated wind,  W , .   Phillips  (1956)  suggested  that  for 
simplicity  in  stress  calculations,  the  anemometer  wind  speed 
can  be  assumed  to  be  70  per  cent  of  the  geostrophic  wind 
speed  \V ,  .   The  stress  was  therefore  taken  to  be 


U    =  P(CD)4'£2KI  W< 


(2.4) 


where  A    =  0 .7 


A  method  for  estimating  C   used  by  Manabe  (1969a)  is 

D 

o 

where   k   =  Von  Karman's  constant, 
o 


(2.5) 


h  =  height  at  which  the  stress  is  to  be  evaluated 
and     z   =  roughness  parameter. 

The  quantity  O   was  assumed  to  be  a  constant  equal  to  4.5 

-3 
x  10   .   This  is  larger  than  values  normally  used  for  this 

coefficient,  but  since  C   depends  on  where  in  the  surface 
layer  the  stress  is  evaluated,  the  appropriateness  of  the 
value  can  best  be  examined  from  the  resulting  momentum  ex- 
change between  the  atmosphere  and  surface.   Values  of  para- 
meters at  level  four  were  obtained  by  linear  extrapolation 
in  p.   The  quantity  \V,  was  thereby  obtained  from 

(2.6) 


W.  =  |  \V,  -  ~W,  . 
4    2   3    2   1 


Smagorinsky  (1963)  derived  a  formula  for  calculating 
the  internal  stress  at  level  two,  by  assuming  that  the 
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stress  was  proportional  to  the  vertical  wind  shear.   A 
slightly  modified  form  of  this  formula  was  used  in  this 
work  and  can  be  written 


*2  =  (pk)2  ^2  (W3  "V* 


(2.7) 


where  (pk)_  is  an  exchange  coefficient  probably  dependent 
on  the  local  Richardson  number.   An  estimate  of  this  from 
Rossby  and  Montgomery  (1935)  and  used  by  Smagorinsky  was 
50  gm  cm   sec   .   Other  values  of  this  coefficient  are 
225  gm  cm   sec    (Palmen,  1955)  and  500  gm  cm   sec   (Riehl, 
1951).   The  choice  of  50  gm  cm   sec    for  this  model  was 
based  on  many  preliminary  numerical  experiments  with  the 
model.   The  quantity  (•"-)  9  is  the  inverse  thickness  of  the 
750-250  mb  layer  and  it  is  equal  to  0.1266  km   .   If  typi- 
cal atmospheric  wind  values  are  used,  the  ratio  of  tT,  to  (C  „ 
is  on  the  order  of  ten  to  one. 

The  horizontal  domain  of  the  model  is  shown  in  Figure  2. 
The  east-west  extent  is  4800  km  and  the  north-south  distance 
is  6000  km.   The  width  was  chosen  as  4800  km  to  coincide 
with  the  wavelength  of  the  x-wave  that  will  be  discussed 
shortly.   The  southern  boundary  is  at  18N  and  the  northern 
boundary  is  at  72N.   The  use  of  the  3~plane  over  such  a 
latitudinal  extent  may  affect  the  results,  but  only  in  a 
quantitative  manner. 

If  the  hydrostatic  approximation,  the  equation  of  state, 

and  the  geostrophic  condition  (ip  =  % —  z)  were  applied  to 

o 
the  first  law  of  thermodynamics,  a  quasi-geos trophic 


23 


6000 


72°N 


18°N 


Figure  2.   Horizontal  domain  of  the 
model.   Distance  in  kilometers. 


version  of  that  equation  would  be  written 


|_(|4,  +(kxV*  •  v<|4)  +  §- 


w  =  -  £-  Q  , 
f  p 
o         o 


(2.8) 


where 


a  = 


R  T    dl  o      o     o 

s  r   s    g  ,         2    -2   -2 

—y L^7 —  +  "p- J  =  2.3  m  sec   cb   , 

P  g  P 


K  -  R/Cp  , 

Q  =  heating  rate  per  unit  mass 


(2.9) 


and 


T  -  temperature  from  the  standard  atmosphere. 
s      r 
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The  above  constant  value  of  0,  the  static  stability,  was 
used  throughout  the  experiment. 

Equation  (2.1)  was  applied  at  levels  one  and  three  and 
Equation  (2.8)  was  applied  at  level  two.   The  following 
equations  were  obtained: 


fj  (v2V  +  Ik  x  v^  •  v(v  tl)  +  BD  ^  - .-gf 


»  -  g 


3t     3t 
i  _ 


3p   3x     3y   1 


a     2                      2  9lN    f  w? 

fj  (V  #3)  +(k  x  V*3  •  V(V  *3)  +  B„  ^  +  -f^ 

n.     3x  3t 

8  3p  L3x  3y  J3 


(2.10) 


(2.11) 


and 

3  A"  V  ,  „     „  »l+  A  „  A-  A   °W2 
( — r— )+  Ik  x  V( = )-V( r—-)-   ~T~ 


3tv  A  p 


A  p 


f  P 


o^2 


Q2  (2.12) 


where   A  p  =  500  mb . 


Following  Lorenz  (1960),  lj;   and  ip_  were  defined  as  follows 


* 


*1+  ^3 


M 


and 


(2.13) 


V  *3 
*M  =  —2 


Then  Equations  (2.10)  and  (2.11)  were  added  and  divided  by 
two  to  yield: 


-^(V2^)    +»k   x    V^M    -V(V2^M)    +  Ik    x    V^T    •    V(V2^T) 


+    3 


di> 


M 


o    3x 


3t  3t 

A  P     L3x  -r-^ J4    * 

3y 


(2.14) 
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Equation  (2.12)  was  solved  for  co2  .   Equation  (2.11)  was 
subtracted  from  Equation  (2.10)  and  the  result  divided  by 
two.   The  quantity  0)9  was  replaced  and  the  following  was 


ob  tained : 

|^(V2-    y2HT   +lkx    VtjJM    •    V(V2-    y2)^T   +  Ik    x    Vih 


V(V  *M> 


d\p 


+   3 


M 


o  3x 


2 

2f   ^2    Ap 
o 


j2fe 


3x     3t 
X £ 


3y  J2 


3t     8t 


9x 


3y  J4 


(2.15) 


2f 


where   y 


A  pa 


=  3.69  x  10"16  cm"2 


To  reduce  computing  requirements,  the  model  was  further 
simplified  by  expressing  the  dependent  variables,  ty       and 
\\)       as  functions  limited  to  one  wave  in  the  zonal  direction. 
The  wave  chosen  was  of  the  predominate  baroclinic  scale. 
Its  wavelength  was  4800  km,  the  zonal  width  of  the  model 
atmosphere.   Any  atmospheric  model  that  is  developed  with 
the  intention  of  reasonably  portraying  air-sea  interaction 
must  include  baroclinic  disturbances  containing  wave  num- 
bers of  about  five  to  eight.   The  baroclinic  processes  that 
are  responsible  for  the  resupply  of  kinetic  energy  to  the 
atmosphere  are  reasonably  represented  by  the  wave  number 
chosen.   Phillips  (1956)  and  Smagorinsky  (1963)  have  demon- 
strated in  experiments  with  no  forcing  functions  in  the 
east-west  direction  that  much  of  the  disturbance  energy  is 
isolated  in  a  single  predominate  wave  number. 
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By  transforming  the  system  of  equations  to  wave-number 
space  in  the  zonal  direction,  grid  points  were  eliminated 
in  the  x  direction  and  the  requirement  for  relaxation 
in  order  to  solve  for  ip  was  removed.   This  type  of  approach 
had  further  computational  advantages  in  that  no  linear 
truncation  errors  arose  in  x-space.   Also,  the  fact  that 
there  were  no  nonlinear  truncation  errors  in  x-space  re- 
duced the  possibility  of  nonlinear  instability.   Therefore, 
errors  caused  by  aliasing  due  to  x-dif f erencing  were  elim- 
inated.  The  fields  were  defined  as  follows: 


and 


ty      =  E(y,t)  +  A(y,t)cos  kx  +  B(y,t)sin  kx   (2.16) 


\p      =  F(y,t)  +  C(y,t)cos  kx  +  D(y,t)sin  kx   (2.17) 


where  A,  B,  C,  and  D  are  Fourier  amplitudes  of  the  disturb- 
ance, and  E  and  F  are  the  zonal  means.   The  quantity  k  is 
the  zonal  wave  number  and  a  value  of  (2tt/4800  )  km    was 
chosen . 

In  the  initial  development  of  the  atmospheric  model  the 
heating  was  simply  proportional  to  the  air-sea  temperature 

difference.   A  sea  surface  temperature,  T    ,  was  specified 

sea 

from  climatological  data  and  held  constant  in  time.   Later 
in  the  work  a  more  complete  form  of  the  heating  was  used 
and  this  will  be  discussed  in  a  later  section  of  this  thesis 

The  mean  temperature  of  the  layer  one-three  is  propor- 
tional to  ip  .   By  using  the  geostrophic  approximation  and 
the  equation  of  state,  the  following  equation  was  obtained: 
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2f 

T  '  =  ! 

2      R   rT 


'K 


(2.18) 


where  T  '  is  a  departure  from  the  mean  temperature  at  level 
two.  To  obtain  physically  reasonable  values  for  the  atmos- 
pheric temperature  at  level  two  a  constant  was  added  to  T  '  . 


Then 


2f 


R   rT 


iJ/_  +  250.3  . 


(2.19) 


By  using  the  standard  lapse  rate,  T_  was  extrapolated  down- 


ward to  obtain  the  temperature  at  level  four,  T. ,  and 


8T 


T,  =  T 


(A  z) 


(2.20) 


4     2    8Z 

where  A  z  is  the  standard  thickness  of  the  1000-500  mb 
layer  and  dT    /dz    is  the  standard  lapse  rate.   Therefore, 
T,  only  differs  from  T   by  a  constant. 

Since  the  static  stability  was  held  constant,  the  heat- 
ing at  level  two  was  proportional  to  the  air-sea  temperature 
difference  and  was  expressed  as 


.£_ 


H,   (T 


T,) 


2    A  p    1    sea     4 


(2.21) 


where  H   is  a  factor  derived  by  Haney  (1971)  and  is  a  slow- 
ly varying  function  of  latitude.   In  this  work  a  value  of 

-2    -1 
55  cal  cm   deg    was  used  and  held  constant  in  both  latitude 

and  time . 

When  Equations  (2.16)  and  (2.17)  were  substituted  into 

Equations  (2.14)  and  (2.15)  and  terms  of  wave  number  2k  and 

higher  were  neglected,  equations  in  the  Fourier  coefficients 

were  obtained.   The  various  terms  were  separated  and 
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coefficients  of  the  cosine  terms,  sine  terms  and  terms  in- 
dependent of  x  were  equated.   This  yielded  six  predictive 
equations  for  A,  B,  C,  D,  E,  and  F.   The  derivation  of  these 
equations  is  contained  in  Appendix  A. 
The  equations  are: 


3  rd     A 

3tl.  2 
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2         3E  32B    3F  82D 
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B.   METHOD  OF  SOLUTION 

1 .   Finite  Difference  Scheme 

Centered  time  differences  were  used  for  all  terms 
except  those  involving  friction  and  heating.   Those  terms 
were  computed  at  time  (t-  At  )  .   A  forward  time  step  was 
utilized  to  start  the  marching  process.   A  finite  differ- 
ence scheme  developed  by  Matsuno  (1966)  was  used  every  24 
time  steps.   This  was  a  two-step  iteration  to  simulate  the 
backward  difference  method.   This  scheme  was  applied  in 
this  situation  to  eliminate  the  computational  mode  that 
would  arise  from  using  the  centered  time  differencing  ex- 
clusively and  also  helps  eliminate  nonlinear  instability. 
A  time  step  of  30  minutes  was  used. 
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Since  a  wave  number  representation  was  used  in  the 
x  direction,  finite  differencing  was  not  required  in  that 
direction.   However,  centered  space  differences  were  taken 
in  the  y  direction.   A  uniform  grid  distance  of  200  km  was 
used  for  Ay . 

2  .   Boundary  Conditions 

Free-slip  walls  existed  at  the  north  and  south 
boundaries.   The  boundary  conditions  at  these  walls  were 

A,  B,  C,  D  =  0     at     y  =  0,   y  =  w         (2.28a) 
and 


3F(37}  =  °>  TE^     =  0  at  y  =  0,  y  =  w. 


(2.28b) 


Condition  (2.28a)  was  imposed  to  force  the  meridional  flow 

to  equal  zero  and  correspond  to  the  statement  that  the 

normal  geostrophic  velocity  vanished  at  the  wall.   From 

3u 


the  zonally  averaged  equation  of  motion, 


3t 


-S.  = 


0  at  the 


north  and  south  walls.   The  boundary  condition  (2.28b)  en- 
sures that  this  was  true  and  was  also  a  requirement  for 
proper  energy  conservation. 
3 .   Initial  Conditions 

Initial  conditions  of  a  mean  zonal  current  with  a 
superimposed  disturbance  were  specified.   The  vertical 
shear  and  hence  the  thermal  wind  were  near  the  critical 
value  for  baroclinic  instability.   The  initial  disturbance 
had  a  wavelength,  previously  discussed,  and  a  structure 
such  that  when  superimposed  on  a  baroclinic  current  with  no 
horizontal  shear  it  had  the  maximum  growth  rate. 
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and 


The  initial  conditions  were: 

200   .   ,-ny.  2     -1 

A  =  — - —  sin  (—*-)     cm   sec 

f         w  ' 


B  =  C  =  D  =  0  , 

E  =  1200  -  cos  (r^-)     cm2  sec"1 

IT         W 


/TTy  ^    2 


F  =   600  —  cos  ( — *-)  cm   sec 

it       w 


(2.29a) 
It  was  further  desired  that  the  mean  surface  wind  vanish 


initially,  which  can  be  achieved  by  requiring  that 

'  =  "  T-t^vr  -  2i/j  ]  =  -  !_[E-2F]  =  0     (2.29b) 

9y   M      x       dy 


U.  =  - 

4      dy 


or 


E  =  2F. 
4  .   Time  Integration  Process 

Different  methods  for  solving  the  time  tendencies 
of  the  Fourier  coefficients  were  utilized.   A  Gauss-Jordan 
elimination  method  was  used  to  solve  equations  (2.22), 
(2.23),  (2.24),  (2.25),  and  (2.27).   Richtmyer  (1957)  dis- 
cussed the  use  of  such  a  procedure  for  solving  a  set  of 
equations  that  had  all  the  elements  of  its  corresponding 
matrix  vanish  except  those  on  the  three  main  diagonals. 
These  equations  had  that  property.   The  method  of  solution 
for  Equation  (2.27)  was  discussed  by  Phillips  (1956)  and 
modified  to  ensure  that  the  equation  satisfied  the  boundary 
conditions.   Equation  (2.26)  was  solved  by  a  direct  march- 
ing process  utilizing  the  proper  boundary  conditions. 


32 


C.   PRELIMINARY  EXPERIMENTS 
1 .   Hadley  Regime 

To  examine  the  effect  of  the  disturbances  on  the 
model  atmosphere,  an  experiment  was  conducted  which  should 
have  produced  a  single  Hadley  type  cell  in  the  meridional 
circulation.   A  linear  sea  surface  temperature  was  imposed, 
with  warm  temperatures  in  the  south. 

The  amplitudes  of  the  disturbance  were  set  equal  to 
zero  (A,  B,  C,  and  D)  and  Equations  (2.22)  through  (2.27) 


reduced    to : 
„2 


9y 


(|1)     .    1_E_JL_    pc       |vv.|     (2    |I-    M)] 
2     v9t  3y     i  p         D    '     41  9y         9y 


and 

r32      2i  3F     9  r  S  /t    ^\        /P8n   8F 
[~ 2  "  ^  9T  =   97[  A?  4(pk)2  <p  >2  "37 


_£ 


2 
H  Kg 

2f  A  p 
o 


Hl  [Tsea  "  (X  (F  +  Cl>)] 


(2.30) 


(2.31) 


For  simplicity  the  following  definitions  were  made: 


-1     Ap  PCD 

S2  E  iV  A(pk)2  (^2  ■ 


and  S 


_  y  Kg 

2f  A  p 
o 


(2.32) 


An  assumption  that  a  steady  state  condition  existed 

in  the  model  was  made,  i.e.  tt— (  )  =  0,  yielding 

d  t 


0  -  f-is,  |\V.|  (2|£  -  M)]. 
dy   1  '  4 '    dy    dy 


(2.33) 
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Since  S..  was  not  equal  to  zero, 


lw, 


=  0 


and 


3E 

dy 


8y 


(2.34) 


and  Equation  (2.31)  can  be  written 

2f 


82F 

3y 


sqhi  lT 

3  1    sea 


R 


(F  +  C 


1>]  "  ° 


(2.35) 


or 


S2  „  2    S3H1   R 


8y 


2f  2f 

F  -  S_H..T     +  S_H..  — ^  C.  =  0 
3  1  sea     3  1   R    1 


(2.36) 


Since  F  represents  the  mean  temperature  throughout 


the  layer,  the  zonal  mean  thermal  wind  can  be  defined  as 

9y  " 


U, 


(2.37) 


Equation  (2.36)  was  solved  analytically,  using  the  boundary 


condition  that  U   =  0  at  y  =  0,  W. 
The  solution  is 
U 


F  = 


where 

U 
and 


5 

8T 

3y 


1/2 


w. 


sea 


(V  _  ™) 1 

^    2;    (S4S5) 

S„ 


1/2 


H. 


sinh[(S4S5)    (y-y) 
cosh[(SAS5)1/2  |] 


S   "1   * 
b2   X 


2f 


(2.38) 


5  "   R 

This  solution,  along  with  the  condition  that  A  thru 
D  equal  zero  and  E  equals  2F,  was  used  as  an  initial  state 
to  expedite  reaching  a  steady  state.   The  initial  mean 


thermal  wind  obtained  from  (2.38)  was 
U 


U, 


cosh[(S4S5)1/2(y-  f)] 
cosh[(S4S5)1/2  |] 


U 


(2.39) 


where  U   =  0.00433  degrees  per  kilometer  and  the  change  in 
the  thermal  wind  with  time  was  examined. 
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The  complete  equations  were  integrated  in  time 
starting  from  this  initial  state.   After  40  days,  values  of 
U   were  within  two  percent  of  the  initial  value.   The  mean 
meridional  circulation  taken  from  the  numerical  solution  is 
shown  in  Figure  3a.   It  consisted  of  a  single  cell  with 
northward  flow  at  the  upper  level  and  southward  flow  in  the 
lower  level.   The  vertical  motion  field  consisted  of  upward 
motion  in  the  southern  part  and  downward  motion  in  the  north 
The  central  portion  displayed  weak  and  variable  vertical 
motion  patterns  and  were  interpreted  as  being  meaningless. 
It  follows  from  the  equations,  that  the  quantities  A 
through  D  will  remain  zero.   The  method  by  which  vertical 
motions  and  the'mean  meridional  motions  were  calculated  is 
described  in  Appendix  B. 

A  similar  experiment  was  made  with  the  same  linear 
sea  surface  temperature  profile  as  before  but  with  the 
initial  conditions  given  by  Equation  (2.29a).   This  initial 
condition  primarily  differs  from  the  previous  one  in  that 
it  contains  a  finite  amplitude  disturbance  with  wave  number 
k  in  the  barotropic  part  of  the  motion  (A  ^  0) .   This  in- 
tegration required  more  computer  time  to  reach  a  quasi- 
steady  state.   A  three  cell  structure  as  shown  in  Figure  3b 
developed  after  230  days.   The  northward  and  southward 
motions  were  stronger  and  the  vertical  motion  patterns 
were  more  intense  and  well  defined. 
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Figure  3.   a.   Meridional  circulation  arrived  at  in  Hadley 
regime  experiment.   b.   Meridional  circulation  of  atmos- 
phere with  disturbances  included. 


A  quantitative  comparison  was  made  of  velocities 
obtained  in  the  two  experiments.   Table  1  lists  the  values 
of  the  upper  level  meridional  velocity  in  cm  sec    and  the 
mean  zonal  velocities  in  m  sec    at  levels  one  and  three. 
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LATITUDE 

MEAN  MERIDIONAL  VE- 
LOCITY UPPER  LEVEL 

ZONAL 

VELOCITY 

LEVEL  1 

LEVEL 

3 

HAD 

DISTURB . 

HAD 

DIST. 

HAD 

DIST. 

72.0 

0  .0 

0.0 

0.0 

0.0 

0.0 

0.0 

68.4 

2.9 

3.4 

20.5 

6.6 

8.0 

-0.1 

64.8 

3.5 

12.2 

18.9 

6.0 

6.8 

-1.9 

61.2 

3.6 

16  .0 

19  .2 

5.0 

6.8 

-1.1 

57.6 

3.5 

6.0 

19.4 

16.3 

6.7 

5.1 

54.0 

3.6 

-5.3 

19.2 

27  .3 

6.8 

11.7 

50.4 

3.6 

-13.7 

19.3 

34  .9 

6.7 

15  .6 

46.8 

3.6 

-22  .2 

19.3 

36.8 

6.8 

16.6 

43.2 

3.6 

-28.8 

19.3 

33.3 

6.7 

14.6 

39.6 

3.6 

-23.8 

19.2 

24.7 

6.8 

9.8 

36.0 

3.5 

-8.6 

19.4 

14  .3 

6.7 

3.7 

32.4 

3.6 

6.4 

19.2 

7.2 

6.8 

-0.3 

28.8 

3.5 

10.1 

19.3 

7.7 

6.8 

-0.8 

25.2 

3.5 

3.1 

18.9 

13.0 

6.8 

0.7 

21.6 

2.9 

1.4 

20.6 

11.3 

8.0 

1.8 

18.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

Table  I.   Comparison  values  of  meridional  and  zonal  velo 
cities  from  Hadley  experiment  and  disturbed  atmosphere 
experiment.   Positive  values  are  northward  or  eastward 


flow.   Upper  level  is  500-0  mb  layer.   Units  ar 


e  m  sec 


-1 


37 


2 .   Energy  Examination 

Equations  for  the  components  of  potential  and  kine- 
tic energy  and  energy  transformations  in  terms  of  the  Four- 
ier coefficients  are  as  follows:   (The  derivation  of  these 
equations  is  shown  in  Appendix  C.) 


uz     2 
;  Pz  dy  =  2   7  F   dy  » 

/  Pe  dy  =  |  f    (C2  +  D2)  dy  , 

/Kz  dy  -I   /  [(|f)2  +  (ff)2]  dy  , 

r   v      a  1       r    r/-9A<*2    /9B^2    rdc\2         rdJ)\2 

f    Ke  dy  =  j      /  I  ( 37)   +  (37)   +  (37)   +  (37) 

+  k2  (A2  +  B2  +  C2  +  D2)  ]  dy  , 

2f  , 

-^(F  +  C1)]F   dy  , 


C   — /  [|  y2(BC  -  AD)  -£rJ  dy  , 


sea 


82B 


•H.-'I^V 


9y 

-  b^A  +  c 


+  -x — (A -  D 

8y  dy2 


3yz 

82A 


+  C 


9y 


92D 

9y2 

32B 

9y2 


-  D 


32C. 

9y2' 


-  B^)]  (dy, 


9y 


Dei=^  jfe  A(pk)2<£f>2^( 


+  D( 


&.    -    Ck2) 

9y2 

d2D 


9y 


-  Die  )]   dy  , 


es 


JS?  PS   [(A"C)  (lvV4|k2(A-2C) 


+  I7[lnl<2l7-|7)])+  (B"C) 


(2.40) 
(2.41) 
(2.42) 

(2.43) 
(2.44) 
(2.45 


(2.46) 


(2.47) 


|\V,|k2(B-2D)    +   i-I  |\V.  I  (2|£   -    l^)] 

1     4  '  dy     '     4 '        dy  dy 


dy,     (2.48) 
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(2.49) 


Dzs  "  -'{fF  PCD  [<E"F>  lTliVAl  ^2ff  -  !f>]]  }^.(2.50) 

(2.51) 


2f 


Cz  Ixa  u,  F]  dy 


and 


2f 

=  ""'         [TT  (wcc   +  %D)]   d*  • 


(2.52) 


The  relationship  of  these  quantities  can  be  expressed  by 
the  following  equations: 


|-  /   Pdy   =G   -C   +  C 
d  t       z         z      p      z 


£-'/   Pdy   =  C   -C 

dt      e        p     e 


~-  /   K  dy   =  C   -C.-D.-D 
dt      e        e     k     ex     es 


and 


(2.53) 
(2.54) 
(2.55) 


yr  /   K  dy 
d  t        Z 


-  C   +  C,   -  D  . 
z     k     zi 


zs 


(2.56) 


If  Equations  (2.53)  to  (2.56)  are  summed,  the  following 
equation  is  obtained: 


~  f       (P+K)  dy=G   -D.-D    -D.-D 
dt  z     ei     es     zi     zs 


(2.57) 


and  shows  that  the  time  rate  of  change  of  total  energy  in 

the  system  is  equal  to  that  generated  minus  that  dissipated 

The  quantity  C   represents  the  conversion  of  mean 
p 

potential  energy  to  eddy  potential  energy  that  is  accom- 
plished by  the  horizontal  transport  of  sensible  heat  by 
the  eddies.   The  temperature  in  the  layer  one-three  is  pro- 

ri  F 

portional  to  \p       and  k(BC  -  AD)  -r—  is  proportional  to  v'T' 
1  d  y 
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9T 

3y 


The  quantity  k(BC  -  AD)  represents  v'T'  which  is  the 

8T 


transfer  of  heat  across  a  latitude  circle,  and  since  ■= — 

dy 

tends  to  be  negative,  k(BC  -  AD)  must  be  positive  to  convert 

P   to  P  . 
z      e 

The  Term  C   represents  the  conversion  of  eddy  poten- 
tial energy  to  eddy  kinetic  energy  by  the  vertical  circula- 
tion in  the  zonal  plane.   The  conversion  depends  on  the 
correlation  between  to'  and  a'  (a'  is  a  departure  from  a 
mean  specific  volume  of  air).   Since  to '  a '  is  proportional 
to  to  '  ip  '  ,  to   and  C,  and  to   and  D  must  be  negatively  corre- 
lated to  convert  P   to  K  . 

e      e 

The  quantity  C,  represents  the  conversion  of  eddy 
kinetic  energy  to  zonal  kinetic  energy  and  is  dependent  on 
the  correlation  between  the  mean  zonal  wind  and  the  meri- 
dional convergence  of  the  eddy  momentum  transport.   The 

j\ 

quantity  u  -r — (u'v')  is  proportional  to 

y         2       2       2       2 

k-r— (A =r  -  B =•  +  C =■  -  D ;r) 

dy  r,    2     .2     -  2       2 

dy      9y      dy      dy 


dy   ~  2 


di!a  .  Bi!c  +  ^^ 

9y      3y      dy 


o-     3E    ,  9F     ,     ,         .     3/t.N       •, 
Since  -r —  and  t: —  tend  to  be  negative,  -r — (u  v  )  must  be 
9y     9y  3y 

positive  in  order  to  have  conversion  of  K   to  K  . 

e      z 

The  term  C   represents  conversion  of  zonal  poten- 
z  r 

tial  energy  to  zonal  kinetic  energy  by  the  mean  meridional 

circulation.   It  depends  on  the  correlation  between  to  and 

a,  which  is  proportional  to  to  ip     .   This  conversion  is  small 

in  nature  (  Oort,  1964)  and  requires  a  positive  correlation 

to  convert  K   to  P  . 
z      z 
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The  quantity  G   is  the  generation  of  zonal  potential 

energy  through  diabatic  process.   In  this  model  the  heating 

is  a  function  of  the  air-sea  temperature  difference.   The 

correlation  of  T  and  Q,  or  F  and  Q,  must  be  positive  to 

generate  zonal  potential  energy.   This  means  that  heating 

of  relatively  warm  air  and  cooling  of  relatively  cool  air 

occurred  in  the  model. 

The  dissipation  terms  D  .  ,  D   ,  D  .  and  D    provide 

eieszi       zs 

an  energy  loss  mechanism  for  eddy  and  zonal  kinetic  energies 
The  second  subscript,  s  or  i,  refers  to  dissipation  by  sur- 
face or  internal  friction  respectively. 

An  interesting  way  to  examine  the  system  energy  is 
by  considering  an  energy  flow  diagram.   Figure  4  depicts 
such  a  diagram.   The  boxes  contain  time  rate  of  change  of 
energy,  and  the  arrows  show  the  direction  of  the  flow  of 
energy  if  the  transformation  is  positive. 

The  model  equations  were  integrated  until  a  statis- 
tically steady  state  was  reached.   Energy  calculations  were 
made  and  instantaneous  values  for  the  quantities  just  dis- 
cussed are  shown  in  Figure  5 . 

The  diabatic  heating  from  the  sea  surface  creates 

-2 
zonal  potential  energy  that  amounts  to  1.27  watt  m   .   The 

conversion  from  zonal  potential  energy  into  its  eddy 

-2 
counterpart  occurs  at  the  rate  of  1.32  watt  m   .   The  eddy 

potential  energy  is  converted  to  eddy  kinetic  energy  through 

-2 
the  eddy  circulation  at  a  rate  of  1.34  watt  m   .   The 
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Figure  4.   Energy  flow  diagram.   Flow  is  in  the  direction 
of  the  arrows  if  the  transformation  is  positive. 


conversion  of  eddy  kinetic  to  zonal  kinetic  energy,  which 

maintains  the  zonal  currents  against  frictional  dissipation 

-2 
is  0.63  watt  m   .   The  dissipation  of  eddy  kinetic  energy 

-2 
by  surface  friction  is  0.79  watt  m   .   This  is  substantially 

larger  than  dissipation  by  internal  friction  and  the  dissi- 
pation of  zonal  kinetic  energy  by  surface  friction. 

It  should  be  noted  that  the  change  of  potential 
energy,  both  zonal  and  eddy,  was  very  small  while  the  change 
in  kinetic  energy  was  substantial.   The  atmospheric  model 
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Figure  5.  Flow  diagram  of  the  model  atmosphere  energy 
Energy  and  energy  transformation  units  are  in  watt  m~2 
(103  erg  cm-1  sec-1). 


had  reached  a  statistically  steady  state  at  the  time  the 
energy  computations  were  made,  but  the  energy  continued  to 
fluctuate  about  a  mean  value.   This  undoubtedly  influenced 
these  calculations. 

An  exact  balance  was  not  achieved  throughout  the 
entire  system.   It  is  believed  that  these  very  small  dis- 
crepancies were  due  to  truncation  or  roundoff  errors  and 
should  not  be  considered  a  serious  defect  in  the  calcula- 
tions.  Figure  6  is  included  to  show  that  these  values  are 
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Figure  6.   Flow  diagram  of  atmospheric  energy  as  con- 
structed by  Oort(1964).   Energy  units  are  in  105  joule 
m~2  ;  energy  transformation  units  are  in  watt  m~  . 


not  only  qualitatively  correct  but  also  quantitatively 
acceptable.   Figure  6  shows  the  atmospheric  energy  flow 
diagram  as  constructed  by  Oort  (1964).   The  values  en- 
closed in  the  boxes  are  no  longer  time  rates  of  change  of 
energy  but  the  amount  of  energy  present. 

It  is  satisfying  to  note  that  the  transformations 
in  the  model  atmosphere  compare  favorably  with  those  of 

Oort.   The  conversion  C   and  C   agree  very  well  while  C 

p       z  b  J  e 

and  C   lie  just  outside  the  extreme  values  as  proposed  by 
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Oort.   The  dissipation  of  eddy  kinetic  energy  is  much  small- 
er than  Oort's.   Since  his  value  itself  is  considered  too 
small,  this  points  out  a  weakness  in  the  dissipative  mechan- 
ism of  the  model.   This  may  be  due  to  the  lack  of  horizontal 
diffusion  terms  in  the  model  equations.   It  was  decided  that 
diffusion  was  unnecessary  in  the  model  and  this  decision  may 
be  reflected  in  the  energy  dissipation.   Also,  energy  loss 
through  nonlinear  interaction  between  waves  has  been  elimi- 
nated in  this  model  by  the  single  wave  formulation. 
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III.   THE  OCEAN  MODEL 

A.   DESCRIPTION  AND  EQUATIONS  OF  THE  MODEL 

The  ocean  model  is  similar  to  that  studied  by  Bryan 
(1969)  and  Haney  (1971a).   It  is  a  baroclinic  model  based 
on  the  primitive  equations  in  which  steady  state  solutions 
are  sought  by  direct  numerical  integrations  of  an  initial 
value  problem.   The  circulation  produced  in  the  ocean-only 
model  is  primarily  driven  by  surface  forces  that  simulate 
wind  stresses,  and  horizontal  pressure  gradient  forces  due 
to  internal  stratification.   The  horizontal  stress  exerted 
on  the  ocean  surface  by  large  scale  wind  systems  causes  the 
surface  water  to  move  as  an  Ekman  layer  and  drift  to  the 
right  of  the  stress  in  the  Northern  Hemisphere.   This  causes 
the  curl  of  the  wind  stress  to  become  important  as  it  causes 
convergence  or  divergence  of  the  Ekman  layer  transport  and 
corresponding  vertical  velocities  below  the  Ekman  layer. 

Heating  in  low  latitudes  and  cooling  in  high  latitudes 
causes  a  thermal  circulation  to  develop.   A  thermal  current 
is  set  up  where  there  are  easterlies  in  the  south  and  west- 
erlies in  the  north.   These  currents  converging  into  the 
continental  boundaries  contribute  to  the  resulting  circula- 
tion . 

The  rectangular  ocean  basin  is  everywhere  1700  meters 
deep  and  is  4800  km  wide  and  6000  km  in  length.   To  coin- 
cide with  the  atmospheric  domain  the  basin  extends  from  18 
degrees  to  72  degrees  north  (Figure  2) .   As  in  the  model 
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atmosphere,  the  (3-plane  approximation  is  used.   Veronis 
(1963a, b)  has  shown  that  the  p-plane  approximation  is  satis- 
factory for  the  intermediate  scale  motions  in  middle  lati- 
tudes.  This  formulation  is  also  used  to  maintain  consistency 
with  the  atmosphere  and  prepare  for  future  joining  of  the 
two  models.   The  3~plane  approximation  in  both  models  elimi- 
nates the  need  for  mapping  and  reduces  the  possibility  of 
confusion  that  might  arise  if  an  ocean  model  in  spherical 
coordinates  were  joined  with  an  atmospheric  model  using  car- 
tesian geometry. 
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Figure  7.   Vertical  structure  of  the  ocean  model. 
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The  vertical  structure  of  the  three-level  model,  shown 
in  Figure  7,  is  changed  considerably  from  that  described  by 
Haney.   The  three  levels  are  placed  at  100,  500,  and  1300 
meters  below  an  undisturbed  surface.   Temperature  and  hori- 
zontal velocities  are  predicted  at  these  levels  while  the 
vertical  velocity  is  diagnos t ically  computed  at  the  inter- 
mediate levels.   The  vertical  layering  is  defined  in  such  a 
way  that  there  is  a  finer  resolution  near  the  surface  than 
near  the  bottom. 

The  Boussinesq  approximation  (variations  of  density  are 
neglected  except  in  terms  involving  gravitational  accelera- 
tion) and  the  hydrostatic  approximation  are  applied.   In 
the  equation  of  state,  the  dependence  of  density  on  salinity 
is  neglected.   Density  is  assumed  to  be  a  linear  function 
of  temperature  only. 

The  governing  equations  are 
du      1   3d     .  .  1 


lH  =  _  1-  |£  +  vf  +  ^_  F 

dt       p   dx  p    X 


(3.1) 


dv 
dt 


1  12.  f   4.  1 

=  -  -r-*-  -  uf  +  


P0  3y 


Po   y 


(3.2) 


dz 


-    Pg 


vvv  +  |2   =  0 

dz 


(3.3) 
(3.4) 


dT    Q_ 


dt    p  C 


(3.5) 


and 


p  =  po[l  -  a(T  -  Tq)] . 


(3.6) 
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where  F   =  eastward  frictional  force  per  unit  volume, 
x  r  ' 

F   =  northward  frictional  force  per  unit  volume, 

and    Q   =  heating  rate  per  unit  volume. 

The  quantities  F  ,  F   and  Q  represent  the  assumed  con- 

x    y 

tributions  of  sub-grid  scale  motions  to  the  exchange  of 
momentum  and  heat  in  the  ocean.   The  internal  redistribu- 
tion is  accomplished  by  utilizing  a  linear  lateral  diffusion 
of  heat  and  momentum.   Similarly,  the  vertical  mixing  of 
these  two  quantities  by  sub-grid  scale  eddies  is  specified 
by  a  linear  vertical  diffusion  coefficient.   The  vertical 
mixing  of  temperature  is  enhanced  by  the  use  of  an  instan- 
taneous adjustment  mechanism  to  eliminate  unstable  lapse 
rates  . 

For  the  initial  experiments  with  the  ocean  model  a  con- 
stant wind  stress  was  imposed.   Figure  8  shows  the  mean 
zonal  stress  that  was  used  in  this  phase.   The  data  was 
taken  from  Hellerman  (1967). 

Following  Haney  (1971a)  the  downward  heat  flux  at  the 
ocean  surface  is  expressed  as 


K-|^)   1   =  [H   +  HV  (T   -  T    )]/p  C  . 

dz  z=t     lo     la     sea     o 


(3.7) 


The  quantity  H   contains  the  net  downward  flux  of  solar 
o 

radiation  across  the  ocean  surface,  minus  the  upward  flux 
of  longwave  radiation  and  latent  heat  from  an  ocean  surface 
at  a  temperature  equal  to  T  ,  the  prescribed  atmospheric 
equilibrium  temperature.   The  second  term  contains  the  up- 
ward longwave  radiation  and  sensible  and  latent  heat.   The 
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coefficient  H   was  found  by  Haney  to  vary  with  latitude  by 
about  20  percent.   It  was  set  equal  to  a  constant  in  this 
work . 


re  s  s 


Figure  8.   Mean  zonal  stress  used  in  the  ocean  model 
experiments.   The  data  is  from  Hellerman  (1967);  units 
in  dynes  cm~2  < 


B.   COMPUTATIONAL  PROCEDURES 

Since  the  ocean  bottom  is  flat,  w  =  0  is  the  boundary 
condition  on  the  vertical  velocity  there.   Following  the 
method  used  by  Bryan  and  Cox  (1967),  Haney  (1971a)  and 
others,  the  ocean  surface  was  considered  a  balanced  surface 
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at  which  w  =  0.   The  height  of  the  surface  was  not  obtained 
from  the  continuity  equation  so  external  gravity  waves  were 
eliminated,  which  permitted  a  much  larger  time  step  to  be 
used.   However,  the  technique  also  removed  the  vertical 
mean  divergence  from  all  scales  of  motion,  including  Rossby 
waves,  and  introduced  an  artificial  constraint  on  the  phase 
speed  of  these  waves  (Gates  1968).   The  necessity  for  the 
longer  time  step  to  accomplish  long-term  integrations  made 
this  technique  mandatory.   This  filtering  technique  does 
not  affect  the  thermally  driven  part  of  the  motion. 

A  new  set  of  equations  was  formed  by  subtracting  the 
vertically  integrated  form  of  Equations  (3.1)  and  (3.2)  from 
Equations  (3.1)  and  (3.2)  respectively.   After  eliminating 
the  pressure  deviation  in  favor  of  the  temperature  by  the 
hydrostatic  equation  and  the  equation  of  state,  equations 
for  the  time  rate  of  change  of  the  vertical  shear  current 
(u'v')  were  obtained.   These  equations  have  the  form 
8 


|t(u')  =  F1(u,v,T) 


(3.8) 


and 


f^Cv1)  =  F2(u,v,T), 


(3.9) 


where 


and 


u'(x,y,z,t)  =  u(x,y,z,t)  -  u(x,y,T) 


v*(x,y,z,t)  =  v(x,y,z,t)  -  v(x,y,T), 


(3.10) 


(3.11) 
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and  u  and  v  are  the  vertical  mean  motions  discussed  shortly. 
Equations  (3.8)  and  (3.9)  were  numerically  solved  to  obtain 
new  values  of  the  vertical  shear  current. 

A  vorticity  equation,  from  which  the  vertical  mean  cur- 
rent was  predicted,  was  formed  by  eliminating  the  pressure 
between  Equations  (3.1)  and  (3.2).   The  result  was  vertical- 
ly averaged  and  had  the  form 

(3.12) 


r  i/    ^dz  =  F3(u>v) 


■H 


where  £  is  the  relative  vorticity  and  H  is  the  depth  of  the 
basin.   Vertical  integration  of  the  continuity  equation 
(3.4)  allowed  the  introduction  of  a  stream  function  for  the 
vertically  integrated  velocity  since  the  vertical  mean 
divergence  was  zero.   The  vertically  averaged  vorticity 
therefore  had  the  form 


\r 


Cdz  =  V2^ 


(3.13) 


-H 


and  Equation  (3.12)  can  be  written 
9 


It  V^  =  VU'V)  ' 


(3.14) 

A  relaxation  technique  was  used  to  solve  Equation  (3.14) 
for  the  time  rate  of  change  of  ip .   The  quantity  \p    was  then 
set  to  zero  on  the  boundaries.   A  time  integration  was  then 
performed  to  obtain  the  new  value  of  the  stream  function 
and  then  the  vertically  averaged  velocity  components  from 


i  r° 
»J     u 


dz 


M    -  - 


dy 


U 


-H 


and 


W- 


dz  = 


m    81 


3x 


=  v  . 


(3.15) 


-H 


52 


The  total  current  was  then  obtained  by  adding  the  verti- 
cal shear  components  predicted  by  Equations  (3.8)  and  (3.9) 
to  the  vertically  averaged  components  predicted  by  (3.14). 

To  obtain  new  values  of  the  temperature,  Equation  (3.5) 
was  written  as 

8T 


at  =  F4   (U'V>T) 


(3.16) 


and  integrated  in  time  to  obtain  T.   The  quantities  p  and  w 
were  then  diagnos tically  calculated  from  (3.6)  and  (3.4) 
respectively . 

The  boundary  conditions  on  momentum  and  heat  at  the 
walls  were  those  of  zero  slip  and  perfect  insulation.   The 
ocean  bottom  was  also  considered  to  be  perfectly  insulated. 
As  a  simplification  of  the  model,  bottom  friction  was  neg- 
lected and  free  slip  was  allowed. 

C.   HEATING  AND  FRICTION 

It  was  assumed  that  all  solar  radiation  penetrating  the 
ocean's  surface  was  absorbed  in  the  top  layer  which  was 
300  meters  thick.   Since  the  side  walls  and  bottom  of  the 
basin  were  insulated,  the  only  source  of  heat  was  at  the 
surface  and  was  calculated  from  the  thermal  boundary  condi- 
tion (Equation  3.7).   A  parameterization  scheme  by  which 
the  internal  redistribution  of  heat  by  sub-grid  scale 
phenomena  was  represented  by  a  simple  linear  eddy  diffusion 
of  heat  and  an  instantaneous  convective  adjustment  mecha- 
nism was  used  in  the  model.  The  right  hand  side  of  Equation 
(3.5)  can  be  written 
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2_  .  .   vaT  +  |_(K1I)  +  Mil 


C  . 


p  C     H 
ro 


dz v  3z'     A  t 


(3.17) 


Bryan  and  Cox  (1967)  have  shown  that  the  lateral  mixing 

coefficient,  A  ,  for  the  apparent  temperature  should  be 
n 

less  than  that  for  momentum  given  by  A  .   This  is  due  to 
the  larger  scales  of  turbulent  motion  in  the  ocean  having  a 
tendency  to  take  place  along  isentropic  surfaces  rather  than 
strictly  horizontal  surfaces.   The  lateral  mixing  of  momen- 
tum and  temperature  in  the  model  is  somewhat  dictated  by 
the  resolution  of  the  numerical  grid.   The  coefficient 
should  be  sufficiently  small  so  as  not  to  obscure  the 
lateral  transfer  of  the  quantities.   A  minimum  value  of 

both  A.,,  and  A,T  is  needed  to  satisfy  the  zero  slip  boundary 
M        H  J  r  J 

conditions,  and  to  simulate  the  horizontal  eddy  transport 
of  heat  and  momentum  by  scales  of  motion  too  small  to  be 
resolved  by  the  numerical  grid.   In  this  phase  of  the  work, 
coefficients  were  held  constant  with 

A   =  .25  x  10   cm   sec 

rl 

and 

Q        O  _  "I 

A..  =  1.0  x  10   cm   sec 
M 

The  vertical  diffusion  coefficient,  K,  differs  in  that 
the  exact  physical  processes  by  which  vertical  exchange 

takes  place  within  the  thermocline  is  unknown.   Following 

2     -1 
Bryan  (1969)  a  value  of  1.5  cm   sec    was  used  throughout 

this  s  tudy . 


54 


The  convective  adjustment  mechanism  (6T)  is  difficult 
to  state  mathematically.   This  device  is  a  means  by  which 
the  temperature  structure  in  the  model  was  forced  to  remain 
stable.   After  the  convective  and  diffusive  heating  terms 
were  added,  the  vertical  temperature  profile  was  examined 
for  unstable  lapse  rates.   If  the  temperature  field  resul- 
ted in  an  unstable  lapse  rate,  the  temperatures  were  in- 
stantaneously adjusted  to  the  vertical  mean  temperature  of 
those  unstable  layers.   When  this  first  adjustment  was  com- 
pleted, the  temperature  field  was  again  examined  to  deter- 
mine if  any  other  layers  had  been  made  unstable.   The 
adjustment  process  was  repeated  until  all  layers  were 
gravitationally  stable. 

The  quantity  K- —  was  evaluated  as  follows: 

dz 

K-|^),_  1  =  [H   +  hV  (T   -  T    )]/p  C 


dz    k=y 


lx  a 

T,  -  T, 


sea 


K8z;k  =  l^-  =  K(-^- 

2      z,  -  z 


-) 


and 


Kll>k=2l  "  ^T~^> 


K~).   ,1  =  0. 
dz  k=3y 


(3.18) 


The  last  equality  imposed  the  condition  that  the  ocean 
floor  was  insulated  since  the  vertical  temperature  gradient 
vanishes . 

The  vertical  eddy  diffusion  of  heat  at  level  one,  for 
example,  was  then  defined  to  be 
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5-  (iM) 


K(T  -  T  ) 

IH  +  H.  (T  -  T    )]/p  C 7-^ ^r- 

o    la    sea     o       (zi~  z0) 


(zl/2  "  Zl|) 


(3.19) 
The  heating  term  due  to  lateral  eddy  diffusion  of  heat  in 


(3.16)  may  be  written  as 


AH    V    T    "    A*[3ll3?    +    37(37)]< 


(3.20) 


£\  rp  £\  rp 

The  terms  -z —  and  n —  were  defined  at  points  one-half  grid 

8x      dy 

distance  to  the  east  and  north,  respectively,  of  points  at 
which  T  was  stored  (see  Figure  9).   The  boundary  conditions 
were 


— )    1   = 


II)    3 
3x  i=  j 


and 


8-T-)  1 

max    z 

1  . 


3T\ 

8x;i=i 


max 


(3.21) 


where  i  =  —  is  a  point  one-half  grid  distance  west  of  the 

western  boundary  of  the  grid(located  at  i=l)  and  i    +  tt 

J  o     v  max   2 

is  a  point  one-half  grid  distance  east  of  the  eastern 

boundary  (located  at  i  =  i    ).   These  conditions  force  -z — 

max  dx 

to  pass  through  zero  at  the  eastern  and  western  boundaries 

thereby  producing  thermal  insulation.   Similar  boundary 

condition  on  -r —  at  the  northern  and  southern  boundaries 
dy 

ensure  that  the  normal  temperature  gradient  vanished  at 

the  walls,  making  them  perfect  insulators. 

F       F 

The  frictional  terms  and  —2-      in  Equations  (3.1)  and 

P       P 

(3.2)  are  composed  of  two  parts.   These  are  the  sub-grid 
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scale  vertical  eddy  stress  and  the  sub-grid  scale  lateral 
diffusion  of  momentum.   These  can  be  expressed  as 


and 


F 


3t 

t~  +  AM  V2u 
3z      M 


di 


dz 


2-  +  A„  V2v 


M 


(3.22) 


3  T 

As  with  the  term  K-r — ,  the  quantities  x   and  T   were  defined 

dz  x       y 

at  the  half  integer  vertical  levels .  The  vertical  eddy 

stress  was  prescribed  at  the  surface  and  defined  to  be 
zero  at  the  bottom.   The  zonal  stress  was  defined  at  each 
level  as 


),   1   =  prescribed  surface  stress, 

"o     2 


Tx  (ul  "  U2} 

Po  k=12        (Z1  ~  Z2) 


\  <u2  "  u3} 

Po)k=4   =  K  (22  "  Z3) 


x. 


P  W   "  ° 
ro     2 


(3.23) 


and  the  meridional  stress  was  defined  as 


^o     2 


=  0 


-*)    1 


=  K 


(y-L  -   v2) 

^Zl  "  Z2^ 


T 

~Z)     1 

o   •  2 


=  K 


(v2  '  V3} 
(Z2  "  Z3} 
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T 

ro     2 


(3.24) 


Thus  the  prescribed  stress  was  purely  zonal. 

The  x  component  of  the  term  involving  the  lateral  eddy 
diffusion  of  momentum  can  be  expressed  as 


M        M   9x   dx     3y   3y 


(3.25) 


The  quantities  -r —  and  -^ —  were  defined  one-half  grid  distance 

dx  dy 

to  the  west  and  south  respectively  of  the  velocity  storage 

points.   The  boundary  conditions  were 
3u. 


:), 


3x'i=l    A  x 


3_Us 
8x;i=i 


2u)i=l± 
2 


u  ■ 


max 


A  x  i  =  i 


_  ±_        j 


max   2 


and 


5u> 
8y;j=l 

8u. 
3y;j=j 


max 


IE)    ! 

2_u. 
A  y'j-j. 


(3.26) 


max   2 

These  conditions  impose  the  zero  slip  boundary  condition 
on  the  velocity  at  the  north  and  south  walls  and  the  no 
flux  condition  at  the  east  and  west  walls.   Similar  defini- 
tions for  the  meridional  components  impose  the  zero  slip 
condition  on  the  east  and  west  walls  and  the  no  flux  condi- 
tion on  the  north  and  south  walls. 


D.   FINITE  DIFFERENCE  TECHNIQUES 

A  space  differencing  scheme  that  utilized  a  centered 
differencing  with  a  staggered  grid  was  used  in  the  integra- 
tion of  the  primitive  equations.   When  a  centered  space 


scheme  is  used  with  an  unstaggered  grid,  a  computational 
mode  in  space  arises  that  is  analogous  to  the  computational 
time  mode  that  occurs  with  a  centered  time  differencing 
(leap  frog)  scheme.   The  use  of  the  staggered  grid  not  only 
prevents  the  appearance  of  the  computational  mode  in  space 
but  is  responsible  for  considerable  savings  in  computing 
time  b.ecause  variables  are  retained  at  alternate  grid 
points.   Arakawa  (1966)  has  shown  that  maintaining  integral 
constraints  on  certain  quadratic  quantities  of  physical  im- 
portance is  sufficient  to  guarantee  nonlinear  computational 
stability.   The  space  differencing  used  in  this  model  is  a 
simplified,  energy  conserving  version  of  the  scheme  used  in 
the  UCLA  general  circulation  model,  as  documented  by  Gates 
et  al  (1971) . 

Total  energy  is  conserved  in  this  model.   The  boundary 
conditions  of  zero  normal  flow  through  the  side  walls,  and 
zero  vertical  velocity  at  the  bottom,  guaranteed  that  the 
following  integral  constraints  were  satisfied  by  the  differ- 
ential equations. 

Kinetic  Energy 


-J    J  P0U  +  v  )dv 

v 

Potential  Energy 

-J   pg(z+H)dV 


/ [-pgw  +  p  (uF   +  vF  )]dv.    (3.27) 
J  °      o    x      y 


f   [pgw  -  g(z+H) 


aj£- 


(3.28) 


v 


The  potential  energy  was  taken  to  be  zero  at  the  bottom 


59 


If  (3.27)  is  added  to  (3.28),  a  conservation  equation 
for  the  total  energy  in  the  basin  is  derived  and  is  given 
by 


v  = 


-  f  [\    PQ(u2  +  v2)  +  pg(z+H)d 

V 

/[p  (uF   +  vF  )  -  g(z+H)^]dv. 


(3.29) 


Finite  difference  forms  of  the  prediction  Equations 
(3.8),  (3.9),  (3.12)  and  (3.16)  are  shown  in  Appendix  D. 
Haney  (1971a)  has  shown  that  these  difference  equations 
conserve  total  energy. 

Figure  9  shows  the  staggered  grid  network.   The  velo- 
city components,  u  and  v,  are  located  one-half  grid  incre- 
ment east  and  one-half  grid  increment  north  of  the  points 
where  corresponding  scalar  variables  T,  p  and  \p    are  re- 
tained.  The  grid  distance  is  200  km  in  both  the  x  and  y 
directions . 


Figure  9.   Location  of  variables  in  the  horizontal 
staggered  grid . 
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In  the  time  differencing  scheme,  various  terms  in  the 
equations  of  motion  were  treated  differently.   The  pressure 
gradient  term  was  evaluated  using  a  centered  difference. 
The  Coriolis  term  was  evaluated  by  a  trapezoidal  implicit 
scheme,  and  the  friction  term  was  evaluated  using  a  forward 
time  step  to  ensure  linear  computational  stability. 

The  trapezoidal  implicit  scheme  for  the  Coriolis  acce- 
leration reduces  to  a  form 


n+1     n-1    „  .   r  •  ir   n+l    /  -\       \  •  c      n~l  t 
u     -  u     =  2At[aif  u     +  (l-a)if  u    ]. 


(3.30) 


It  is  easy  to  show  that  this  scheme  is  linearly  computation- 
ally stable  for  all  A  t  if  —  \£  a  £  1  .   In  this  model, 
a  =  0.55  . 

To  supress  computational  modes  in  time  that  can   arise 
when  a  leap-frog  scheme  is  used  exclusively  in  an  oscillat- 
ing system,  a  Matsuno  (Euler-Backward)  step  was  used 
periodically.   This  two-step  scheme  damps  the  computational 
mode  and  permits  the  long  term  integration. 

The  leap-frog  scheme  is  conditionally  stable  and  the 
maximum  time  step  was  determined  by  the  highest  frequency 
governed  by  this  scheme.   External  gravity  waves  have  been 
removed  from  the  model  which  leaves  internal  gravity  waves 
as  the  highest  frequency  waves  present.   The  phase  speed  of 
these  waves  is  Vg'h'  where  g'  is  a  modified  gravity  equal 

to  g £.  .   In  this  model  g»  ~  3.60  cm  sec"2,  H'  =  1200 

o 
meters  and  A  y  =  200  km,  and  the  linear  computational  sta- 
bility criterion  is 
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At  v< 


A  y 


. 375  days . 


This  was  the  time  step  used  in  this  phase  of  the  model. 

E.   INITIAL  CONDITIONS 

A  state  of  rest  was  imposed  on  the  ocean  for  the  start 
of  the  first  integrations.   A  stratified  temperature  struc- 
ture was  selected  so  that  a  minimum  amount  of  computing  time 
would  be  required  to  bring  the  ocean  model  to  thermal  equi- 
librium.  Figure  10  shows  the  vertical  temperature  structure 
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Figure  10.   a)  Initial  temperature  field  at  level  one; 
b)  Same  at  level  two;  c)  Same  at  level  three. 
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The  northernmost  temperature  at  each  level  is  A  degrees 
Centigrade  and  the  vertical  structure  is  isothermal  in  that 
column.   The  temperature  at  level  one  increased  from  AC  at 
the  northern  boundary  to  26C  at  the  southern  wall.   At  level 
two  the  temperature  increases  from  AC  to  8.5C  and  at  level 
three  it  is  a  constant  AC  throughout.   The  stability  varies 
from  isothermal  in  the  northernmost  column  to  very  stable 
in  the  southernmost  column.   No  zonal  variations  were  ini- 
tially imposed . 
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IV.   THE  COMBINED  AIR-SEA  MODEL 

A  combined  model  was  developed  by  joining  the  atmosphere 
and  ocean  models.   The  exchange  of  momentum  and  heat  across 
the  air-sea  interface  is  the  mechanism  upon  which  the  cou- 
pling of  the  two  models  was  based.   An  exchange  of  moisture 
was  not  possible  since  the  atmospheric  model  was  "dry",  but 
a  crude  attempt  at  parameterizing  latent  heat  transport  in 
the  atmosphere  was  incorporated  and  will  be  discussed  in 
the  section  on  heat  exchange. 

The  physical  aspects  of  the  combined  models  remained 
the  same  as  they  were  in  the  individual  models.   The  spatial 
grid  in  the  ocean,  and  time  increments  for  each  model,  re- 
mained the  same.   The  time  marching  process  allowed  the 
time  integration  to  remain  synchronized  and  avoided  a  dis- 
torted response  of  the  ocean  circulation  to  atmospheric 
fluctuations.   The  atmospheric  solution  was  evaluated  at 
grid  points  in  the  x-directlon  in  the  atmospheric  model. 
This  was  done  to  construct  a  complete  stress  field  from  the 
atmospheric  Fourier  coefficients  in  order  to  calculate  the 
air-sea  momentum  exchange.   The  ocean  surface  temperature 
was  transformed  to  wave  number  space  to  calculate  the  air- 
sea  heat  exchange.   The  overall  techniques  of  integration 
used  in  each  separate  model  remained  unchanged  in  the  com- 
bined model . 
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A.   TIME  SYNCHRONIZATION 

One  of  the  problems  that  arose  in  early  air-sea  models, 
such  as  that  developed  by  Manabe  (1969a, b)  and  Bryan  (1969) 
was  the  imprac t icali ty  of  using  a  marching  process  that 
would  keep  the  ocean  and  atmosphere  synchronized.   Primitive 
equation  atmospheric  models  require  time  steps  on  the  order 
of  minutes  and  cannot  be  realistically  integrated  in  sequence 
with  ocean  models  that  have  time  steps  on  the  order  of  hours, 
because  the  computing  time  requirements  are  prohibitive. 
Techniques  have  been  developed  to  simulate  time  synchroniza- 
tion but  little  success  has  been  achieved  in  long-term  one- 
to-one  integrations. 

Since  the  present  atmospheric  model  is  quas i-geos trophic , 
a  time  step  of  one-half  hour  was  utilized.   The  ocean  time 
step  was  nine  hours  as  discussed  previously.   The  resulting 
ratio  of  atmospheric  time  steps  to  ocean  time  steps  was 
eighteen  to  one.   The  marching  process  consisted  of  18 
thirty-minute  time  steps  in  the  atmosphere  followed  by  one 
nine  hour  time  step  in  the  ocean. 

Conservation  of  momentum  and  heat  in  the  exchange  pro- 
cess during  the  integration  was  crucial  to  the  development 
of  a  realistic  interaction  model.   Therefore,  the  quantities 
such  as  the  total  surface  stress  and  heating  calculated  dur- 
ing the  18  atmospheric  time  steps  had  to  equal  the  quantity 
as  it  was  applied  during  the  one  ocean  time  step.   The 
method  by  which  this  was  accomplished  will  be  discussed  in 
the  next  two  sections. 
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B.   MOMENTUM  EXCHANGE 

The  surface  stress  discussed  earlier  and  computed  from 
Equation  (2.3)  was  used  as  an  upper  boundary  condition  for 
the  ocean.   The  calculation  of  this  stress  was  based  on  the 
assumption  that  there  was  no  motion  in  the  ocean.   That 
assumption  is  also  used  in  the  stress  formulation  in  the 
combined  model.   Since  the  velocity  of  ocean  currents  is 
much  smaller  than  the  velocity  of  air  flow,  the  effect  of 
ocean  currents  on  momentum  exchange  can  be  neglected. 
Manabe  (1969b)  used  the  same  argument  in  justifying  his  use 
of  a  similar  formulation  of  the  stress.   Simplicity  and 
computing  requirements  also  entered  into  the  decision  to 
use  this  formulation  here. 

Starting  from  a  prescribed  initial  state,  the  atmos- 
pheric equations  were  integrated  over  18  time  steps.   The 
surface  stress  was  calculated  each  time  and  an  average  was 
computed  to  be  used  for  the  ocean  integration.   A  general 
numerical  procedure  for  carrying  out  the  computation  of 
the  time  mean  surface  stress  is  given  by 

NM 

STXXij  "  55 WV1^  (4-x) 

and  NM 

STXYii  =  Wi[L    <Vijn] '  (4-2) 

J       n  =  l   7 

where  NM  =  ( A  t)      /(At) 

ocean       atmosphere 

and  i  and  j  are  horizontal  grid  indices  having  their  usual 

direct  ions . 
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Since  the  surface  stress  components  in  the  atmospheric 

model,  T   and  T  ,  were  not  calculated  directly,  terms  had 
'   x      y 

to  be  collected  and  values  for  the  stresses  recovered  at 
each  grid  point.   The  stresses  were  computed  in  three  equa- 
tions, (2.22),  (2.23)  and  (2.26). 

The  following  definitions  were  made: 

T  )  .  ,  T  r)  .  =  Contributions  of  the  A  disturbance  to  the 
x  A    Y  A 

zonal  and  meridional  components  of  stress 

computed  in  the  A  equation. 

T  )  ,  T  )   -  Contributions  of  the  B  disturbance  to  the 
x  B    Y  B 

zonal  and  meridional  components  of  stress 

computed  in  the  B  equation. 

T  )  ,  T  )   =  Contributions  of  the  mean  field  to  the 
x  E    YE 

zonal  and  meridional  components  of  stress 

computed  in  the  E  equation. 
To  reconstruct  zonal  and  meridional  stress  fields  from 
these  quantities,  the  terms  must  be  multiplied  by  trigono- 
metric functions  as  follows: 

(t  )..  =  T  )_  +  T  ) .  cos  kx .  +  T  )_  sin  kx .  (4.3) 

xij     xE     x  A        l     xB        l 

and 

(x  )..  -    T  )„  +  T  )  .  cos  kx.  +  T  ) _  sin  kx .  .  (4. 4) 

y  i]     yE     yA        r     yB        l 

The  locations  at  which  the  atmospheric  surface  stresses 
were  calculated  coincide  with  the  points  in  the  ocean  at 
which  the  temperature  and  stream  function  were  calculated. 
Since  u  and  v  were  calculated  at  the  staggered  or  half-grid 
points  in  the  ocean  system  (see  Figure  11),  it  was  necessary 
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to  average  the  stresses  calculated  in  (A. 3)  and  (4.4)  before 
applying  the  stress  to  the  ocean.   The  result  is 


u      1 

x  i,j+l      x  i+l,j+l 


<T*>ij      =  t^Vij    +    <Vi+lfJ    +    co,    ^    +    co,^    _] 


(4.5) 
and 

(ViJ     =  T"V±J  +   (Vi+i.J   +   (Vi,J+i  +   Vi+i.J+i1- 

(4.6) 
STXX  and  STXY  were  then  calculated  from  (t  )U  and  (t  )U 

x  ij      y  ij 

for  application  to  the  ocean.   This  procedure  of  averaging 
over  the  number  of  atmospheric  integrations,  and  applying 
the  stress  at  the  proper  position  was  continued  throughout 
the  duration  of  the  long  term  integrations. 

C.   HEAT  EXCHANGE 

To  satisfactorily  conserve  heat  in  the  air-sea  exchange 
process,  a  time  averaging  procedure  similar  to  that  used  in 
the  momentum  exchange  was  utilized.   The  equations  for  the 
heat  flux  in  both  atmosphere  and  ocean  must  guarantee  that 
the  total  heat  lost  by  the  ocean  equals  that  gained  by  the 
atmosphere . 

The  atmospheric  heating  was  previously  discussed  and 
given  by  Equation  (2.21).   In  the  coupled  model,  that  for- 
mulation was  expanded  to  include  the  net  radiative  heating 
of  the  atmosphere  and  a  crude  specification  of  the  heat 
released  through  precipitation  processes.   The  heat  flux 
into  the  atmosphere  was  given  by 

QA=R+L.P+H1(T     -T.)  (4.7) 

A  1   sea     4 
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Figure  11.   Grid  alignment  of  atmosphere  and  ocean  to  show 
non-coincidence  of  stress  points  in  atmosphere  and  velo- 
city storage  points  in  the  ocean. 
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where  R  =  net  flux  of  radiation  into  the  atmosphere, 
and  L  ♦  P  =  net  flux  of  heat  into  the  atmosphere  due  to 
precipitation.   The  specification  of  R  and  L  •  P  are  des- 
cribed below. 

The  ocean  heat  flux  was  also  modified  from  the  form 
used  in  the  ocean  only  model.   The  necessity  for  modifica- 
tion arose  from  the  fact  that  in  the  ocean-only  model  (like 
the  atmosphere-only  model)  there  was  no  requirement  for  an 
exact  accounting  of  heat  exchange.   The  ocean  could  give 
up  any  amount  of  heat  into  the  non-existent  atmosphere  with 
no  adverse  effects.   Similarly,  when  heat  was  required, 
there  was  an  endless  supply  available  from  the  same  non- 
existent atmosphere.   In  the  combined  model,  the  loss  of 
heat  from  the  ocean  by  evaporation  was  a  specified  function 
of  y  only.   The  quantity  H   in  Equation  (3.7)  was  divided 
into  a  term  representing  the  net  downward  flux  of  short 
and  long  wave  radiation,  S  ,  minus  the  upward  flux  of 
latent  heat,  L»E.   The  equation  for  downward  oceanic  heat 


flux,  Q  ,  was  therefore  written  as 
o 

Q   -  S.T  -  L-E  -  H-  (T     -  T.  ) 

o     N  1   sea     4 

where  S>T   =  net  downward  radiative  heat  flux 
N 


(A. 8) 


and  L"E    E  net  downward  latent  heat  flux. 

Figure  12  gives  a  schematic  representation  of  the  heat  flux 

into  the  combined  model  and  the  calculated  exchange 

Hn (T    -T. ) .   The  arrows  originating  or  terminating  out- 
l   s  e  a   4 

side  the  model  atmosphere  or  ocean  indicate  externally 
specified  heat  fluxes. 
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Figure  12.   Schematic  representation  of  simplified  heat 
flux  into  combined  model.   See  text  for  definition  of 
symbols . 


Annual  mean  values  were  used  in  the  radiation  terms  in 
both  Equations  (4.7)  and  (4.8).   Wetherald  and  Manabe 
(1972)  recently  reported  on  the  effects  of  seasonal  varia- 
tion of  the  solar  radiation  on  the  joint  ocean-atmosphere 
model  developed  by  Manabe  and  Bryan.   Their  model  atmos- 
phere showed  a  marked  warming  in  high  latitude.   Although 
a  similar  result  might  be  expected  in  this  model,  varying 
the  solar  radiation  should  be  the  subject  of  a  separate 
investigation  and  will  not  be  included  in  this  work. 
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From  Equations  (4.7)  and  (4.8),  the  total  heat  flux  into 


the  ocean-atmosphere  system,  Q   +  Q  ,  is  given  by 

a     o 

^  a     o     N 


(4.9) 


When  Equation  (4.9)  was  integrated  over  the  horizontal  do- 
main of  the  model,  the  total  prescribed  heat  flux  into  the 
combined  system  was  obtained.   Time  integration  of  the 
model  cannot  approach  a  steady  state  unless  this  total  pre- 
scribed heat  flux  into  the  system  is  zero.   Also,  the  atmos- 
pheric circulation  is  insensitive  to  the  value  of  the  total 
heat  flux  because  only  an  input  of  available  potential 
energy  impulses  the  circulation.   Since  each  of  the  heating 
components  on  the  right  hand  side  of  (4.9)  represent  essen- 
tially different  physical  effects,  each  term  was  prescribed 
to  give  a  zero  net  heat  flux  when  integrated  over  the  hori- 
zontal domain.   The  meridional  mean  was  removed  from  R  and 

S„  and  the  quantities  L»P  and  L*E  were  made  sinusoidal 
N 

about  zero  in  the  y  direction.   Observed  data  shows  that 

the  ocean  receives  radiation  and  the  atmosphere  loses  it 

at  all  latitudes.   Removing  the  mean  and  using  the  new 

quantities  R'  and  S  '  simulates  the  role  that  the  media 

n 

play  as  both  a  source  and  sink  of  radiative  heating.   Re- 
moving the  mean  from  both  L*P  and  L'E  was  necessary  to 
insure  that  the  net  prescribed  heating  in  both  the  atmos- 
phere (given  by  (4.7))  and  ocean  (given  by  (4.8))  was 
zero.   Thus,  there  is  no  net  prescribed  heat  exchange  be- 
tween the  two  media  -  all  of  the  heat  exchange  is  calcu- 
lated internally  from  the  last  term  in  (4.7)  and  (4.8)  and 
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depends  on  the  air-sea  temperature  difference.   Denoting 
new  quantities  that  have  a  zero  meridional  average  by 
primes,  Equations  (4.7)  and  (4.8)  become 


QA  =  R1  +  L-P1 +  H-  (T     -  T. ) 
^A  1    sea     4 


(4.10) 


and 


Q   =  S  '  -  L-E'  -  H,  (T     -  T. ) . 
^o     n  1    sea     4 


(4.11) 

These  terms  will  be  discussed  in  the  next  sections. 

1 .   Atmospheric  Radiation,  R',  and  Solar  Radiation,  S 
The  observed  annual  mean  net  atmospheric  radiative 
cooling  plotted  as  a  function  of  latitude  is  shown  by  the 
dashed  line  in  Figure  13.   The  data  is  from  London  (1957). 
It  shows  a  meridional  average  loss  of  heat  due  to  radiative 
effects  of  184  langleys  day   .   The  do t-and-dashed  line  in 
this  figure  portrays  the  deviations  from  the  mean  radiation 
As  can  be  readily  seen,  the  deviations  from  the  mean  are 
very  small  and  are  in  all  cases  less  than  10  langleys  day 
These  deviations  were  considered  too  small  to  be  signifi- 
cant and  were  neglected.   Therefore,  R'  =  0  for  future 
calculations  . 

Figure  13  also  shows  the  observed  annual  mean  solar 
radiation  minus  the  mean  upward  longwave  radiation  at  the 
earth's  surface,  S  .   This  quantity  with  its  mean  removed, 
S  ',  is  shown  by  the  dashed  line  in  Figure  14.   The  values 
of  this  quantity  used  in  the  model  are  shown  by  the  solid 
line  in  this  figure.   A  profile  that  was  both  linear  and 
symmetric  was  chosen  and  held  constant  in  time. 
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Figure  13.   Net  mean  atmospheric  radiation,  R,  (dashed 
line);  deviation  from  mean  atmospheric  radiation,  R', 
(do t-and-da shed  line);  annual  mean  net  radiation  arriving 
at  surface,  S  ,  (solid  line).   Data  from  London  (1957). 
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Figure  14.   Deviation  from  mean  net  radiation,  S 


line)  as  observed  by  London  (1957) 

radiation,  S  '  ,  (solid  line). 
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2 .   Latent  Heat  Terms,  L»E'  and  L » P  ' 

The  parameterization  of  the  latent  heat  processes 
could  be  accomplished  in  one  of  three  ways.   Since  the  at- 
mospheric model  contained  no  moisture,  a  direct  exchange 
of  this  quantity  and  the  associated  releases  of  heat  could 
not  be  considered.   One  possible  formulation  was  based  on 
the  assumption  that  the  air  directly  above  the  sea  surface 
was  saturated  so  no  evaporation  existed  (L*E'=  0).   This 
formulation  coupled  with  the  assumption  that  no  precipita- 
tion existed  thereby  eliminated  the  terms  L*E'  and  L'P1 
from  Equations  (4.10)  and  (A. 11).   An  experimental  short 
term  integration  of  the  model  equations  was  completed  using 
this  formulation.   Results  were  unsatisfactory  in  that  it 
became  evident  in  this  short  test  that  an  active  air-sea 
model  could  not  be  developed  with  this  scheme.   The  merid- 
ional  gradient  of  the  sea  surface  temperature  became  weak 
and  the  atmosphere  consequently  became  baroclinically 
stable.   After  a  short  time  this  formulation  was  discarded. 

A  formulation  of  L'E'  and  L»P'  that  was  similar  to 
the  observed  profiles  (Figure  15)  was  used  in  a  second 
short  term  experiment.   The  result  was  similar  to  the  case 
just  described,  only  the  meridional  temperature  gradient 
was  weakened  at  a  more  rapid  rate.   The  solar  radiation 

absorbed  at  the  ocean  surface,  S  '  was  insufficient  to  main- 

N 

tain  an  adequate  meridional  sea  surface  temperature  gradi- 
ent and  the  addition  of  latent  heating  in  the  north  and 
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cooling  in  the  south  further  reduced  the  atmosphere's 
meridional  temperature  gradient.   The  results  of  this  test 
were  also  unsatisfactory.   One  suggestion  for  the  reason 
that  the  ocean  did  not  maintain  a  strong  temperature  gradi- 
ent even  though  the  applied  solar  radiation  was  close  to 
that  observed  in  the  real  atmosphere  is  that  the  model  can- 
not simulate  the  sharp  separation  of  the  western  boundary 
current  from  the  coast  and  the  resultant  tight  gradient  of 
temperature  in  the  mid-latitude  regions.   This  was  probably 
due  to  the  degree  of  linearity  of  the  solution  in  this  area 
(see  Section  VI  A  1)  and  resulted  in  an  oceanic  and  atmos- 
pheric meridional  temperature  gradient  that  was  too  weak 
for  baroclinic  instability.   The  system  as  a  whole  needs 
the  latent  heat  released  in  the  tropics  of  the  real  atmos- 
phere . 

The  results  of  the  above  two  experiments  indicated 
that  the  ocean-atmosphere  model  required  a  somewhat  differ- 
ent partitioning  of  the  heating  than  is  typically  observed 
in  the  real  ocean-atmosphere  system.   The  heating  formula- 
tion adopted  for  the  model  was  equivalent  to  assuming  that 
evaporation  occurred  in  the  ocean  part  of  the  model  and 
that  the  moisture  immediately  precipitated  out  at  the  same 
latitude  in  the  atmosphere  (L*E'  =  L#PT).   This  formulation 
is  somewhat  unrealistic  since,  over  the  latitudinal  domain 
of  this  model,  the  observed  evaporation  exceeds  precipita- 
tion in  lower  latitudes  while  the  reverse  is  true  in  the 
middle  latitudes.   Figure  15  shows  curves  of  the  latent 
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heat  flux  due  to  evaporation  (solid  line),  and  precipita- 
tion (dashed  line),  from  Budyko  (1962).   The  meridional 
average  values  have  been  removed  from  these  values.   The 
sinusoidal  form  of  both  L'E1  and  L*P'  that  was  used  to  speci- 
fy these  fluxes  in  the  model  is  shown  by  the  do t-and-dashed 
line  in  Figure  15.   While  this  formulation  is  quantitatively 
in  agreement  with  the  evaporation  curve,  it  is  a  poor  ap- 
proximation to  the  precipitation  related  heat  release.   How- 
ever, the  greatest  precipitation  is  in  the  equatorial  region 
lying  south  of  the  wall  at  18N.   If  the  hemispheric  preci- 
pitation were  mapped  into  the  domain  of  the  model,  the  pre- 
scribed latent  heat  flux  due  to  precipitation  would  be  more 
in  agreement  with  observed  values.   This  is  considered  to 
be  a  weakness  in  the  heating  formulation  and  is  an  area 
that  could  be  improved  on  in  future  work  with  this  model. 
Nevertheless,  it  was  believed  that  once  the  heating  gradi- 
ent in  the  atmosphere  was  sufficient  to  maintain  it  in  a 
baroclinically  unstable  state,  the  resulting  air-sea  inter- 
actions would  be  realistic.   That  is,  the  interactions 
depend  on  the  characteristic  atmospheric  motions,  independ- 
ent of  the  external  forcing  that  maintain  those  motions. 

Figure  16  shows  the  net  prescribed  atmospheric  heat 
flux,  R1  +  L*P'  (dashed  line)  and  the  net  prescribed  ocean 
heat  flux,  S  '  -  L'E  (solid  line).   Prior  to  its  final 
acceptance  this  form  was  tested  in  a  short  term  experiment 
and  found  to  satisfactorily  contribute  to  the  conditions 
desired  for  an  active  air-sea  model.   This  fact  coupled 
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Figure  15.   Latent  heat  flux  due  to  evaporation  (solid 
line);  latent  heat  flux  due. to  precipitation  (dashed  line) 
as  observed  by  Budyko  (1962)  .   Model  values  of  these 
quantities  (do t-and-dashed  line).   All  have  the  meridional 
averages  removed. 
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Figure  16.   Net  prescribed  atmosphere  (dashed  line)  and 
ocean  (solid  line)  heat  fluxes. 
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with  the  requirement  for  simplicity  led  to  its  use  in  the 

remainder  of  this  work. 

3.   The  Sensible  Heat  Transfer,  H-i  (T     -'  T.  ) 
z x s  e  a 4— 

The  quantity  Hn  (T     -  T.)  will  be  referred  to  as 
n       J       1   sea     4 

the  sensible  heat.   In  the  derivation  of  H   by  Haney  (1971), 
it  was  shown  that  this  quantity  actually  contained  the  net 
upward  flux  of  longwave  radiation  and  sensible  and  latent 
heat  which  occurs  due  to  the  air -sea  temperature  difference. 
Since  the  parameter  H1  varied  little,  a  constant  value  was  . 
assumed.   As  in  the  ocean-only  case  a  value  of  55  langleys 
day    was  chosen. 

As  shown  above,  all  of  the  heat  exchange  between 
the  atmosphere  and  the  ocean  must  be  accomplished  through 
the  sensible  heat  term.   It  is  in  this  term  that  a  strict 
accounting  of  heat  exchange  had  to  be  done. 

The  sensible  heat  flux  applied  during  an  ocean  time 
step  was  the  average  of  the  heat  flux  calculated  during  the 
18  previous  atmospheric  time  steps.   The  equation  utilized 
was 


NM 


HEAT. .  = 


[E   H-  (T     -  T.  )  .  .  ]  , 
nl    sea     4  in  n   ' 
n  =  l  J 


ij  "sea    ~4'ijnJ  (4.12) 

where  NM  =  (At)  ocean/(A  t)  atmos  =  18.   Because  of  the 
x-wave  formulation  used  in  the  atmospheric  model  the  quan- 
tity (T     -  T,)  had  to  be  calculated  in  several  steps. 
J  sea     4 

By  using  Equations  (2.17),  (2.19)  and  (2.20)  an  equation 


for  T,  could  be  written: 
4 

(T  )  .  .  =  (T.F)  .  +(T,C).  cos  kx.  +(T.D).  sin  kx .    ,    (4.13) 
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where 
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(T.F),  =  -~  F,  +  250.3  -  -^r-    (Az) 


4A/j     R    j 


8z 


and 


2f 


(T,D) 
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(4.14) 
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To  be  consistent  with  the  wave  number  representa- 
tion in  the  atmosphere  and  to  facilitate  coupling  the 
models  through  a  single  wave  number,  a  Fourier  analysis  was 
performed  on  the  sea  surface  temperature.   At  each  latitude, 
the  sea  surface  temperature  was  analyzed  in  x-space  to 
yield  a  mean  zonal  sea  surface  temperature  and  the  Fourier 
coefficients  of  its  first  mode.   This  can  be  symbolized  as 

(T    )..  =  (T    ).  +(T    C).cos  kx.+(T    D).sin  kx..(4.15) 
sea  i]      sea  j     sea   j        l    sea   j        l 

The  various  modes  of  the  sensible  heat  flux  in  the  atmos- 
phere were  calculated  from 

H.  (T     -  T. )  -  H.  [ (T     -  T.F)  . 
1   sea     4      1    sea     4   j 


+  (T    C  -  T. C) .  cos  kx. 
sea      4   j        x 


+  (T    D  -  T.D) .  sin  kx . ] 
sea      4   j         i 


(4.16) 


and  were  applied  in  Equations  (2.24),  (2.25)  and  (2.27) 
during  each  time  step.   The  right  hand  side  of  Equation 
(4.16)  was  averaged  over  the  18  atmospheric  time  steps  for 
incorporation  in  Equation  (.4.12).   This  in  turn  was  applied 
in  the  oceanic  heating  formulation. 

The  model  equations  were  integrated  for  a  short 
time  period  (300  days)  and  zonally  averaged  values  of  the 
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sensible  heat  were  calculated.   The  results  are  shown  in 
Figure  17  by  the  solid  line.   The  dashed  line  shows  the  ob- 
served zonal  mean  sensible  heat  flux  as  presented  by  Budyko 
(1962).   These  were  added  to  the  imposed  values  of  L*P'  to 
obtain  the  zonally  averaged  net  atmospheric  heating.   Values 
of  the  observed  sensible  heat  flux  estimated  by  Budyko  were 
added  to  the  observed  values  of  the  radiative  and  latent 
heat  contributions,  and  after  their  meridional  averages 
were  removed,  these  new  quantities  were  plotted  in  Figure 
18.   The  solid  line  shows  the  observed  atmospheric  heating 
over  the  ocean  and  the  dashed  line  shows  observed  atmospher- 
ic heating  over  ocean  and  land  combined.   The  dot-and- 

dashed  line  depicts  corresponding  values  of  L*P!  +  H„ (T 

1   sea 

-  T.)  generated  in  the  model.   In  the  northernmost  segment 
of  the  domain,  the  model  closely  resembles  the  observed 
values  over  the  ocean.   From  a  region  near  45N  to  the 
southern  boundary,  the  model  is  in  closer  agreement  with 
the  data  from  a  land  and  ocean  environment. 

There  are  several  possible  explanations  of  the 
"dip"  in  the  model  atmosphere  heat  flux  near  50N.   The 
formulation  of  the  net  solar  radiation  as  shown  in  Figure 
13  favors  a  sea  surface  temperature  structure  that  would 
be  colder  than  the  atmosphere  in  the  higher  latitudes  and 
warmer  than  the  atmosphere  in  lower  latitudes.   This  would 
lead  to  negative  contributions  from  the  sensible  heat  term 
in  higher  latitudes.   The  curve  dips  sharply  in  the  area 
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Figure  17.  Heat  flux  due  to  sensible  heat  exchange. 
Positive  is  flux  from  the  ocean;  Model  (solid  line); 
Observed    by    Budyko     (dashed    line). 
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Figure  18.   Observed  zonally  averaged  atmospheric  heat  flux 
over  ocean  (solid  line);  observed  zonally  averaged  atmos- 
pheric heat  flux  over  land  and  ocean,  combined  (dashed  line); 
model  zonally  averaged  atmospheric  heat  flux  (do t-and-dashed 
line);  (observed  data  from  Budyko ,  1962). 


where  the  solar  radiation  and  sensible  heat  passes  from 
positive  to  negative  values. 

The  formulation  of  the  latent  heat  released  by  pre- 
cipitation in  the  atmospheric  model  contributes  to  heating 
in  the  south  and  cooling  in  the  north.   The  sharp  dip  in 
the  total  atmospheric  heating  occurs  at  the  latitude  where 
the  prescribed  atmospheric  heating  becomes  negative. 

The  sensible  heat  term  may  contribute  to  this 
abrupt  decrease  in  another  way.   Large  positive  values  of 
this  quantity  are  seen  in  some  regions  and  large  negative 
values  in  others.   The  large  negative  values  are  physically 
unreal  because  the  real  atmosphere  does  not  give  large 
amounts  of  heat  to  the  ocean  through  this  term  over  large 
areas.   These  values  are  considered  departures  from  a  mean 
and  can  be  interpreted  as  small  heating  where  the  values 
are  negative  and  larger  heating  where  they  are  positive. 
The  prescribed  heating  of  the  combined  model  has  already 
been  shown  to  meridionally  integrate  to  zero.   There  can 
be  no  net  heat  flux  by  the  sensible  heat  term  if  a  steady 
state  condition  is  to  be  reached.   If  a  different  formula- 
tion for  the  prescribed  quantities  was  chosen  so  that  the 
net  specified  heat  flux  into  the  system  was  not  zero,  then 

the  net  heat  flux  by  H,  (T     -  T.)  would  have  to  be  adjus- 

J       1    sea     4 

ted.   In  the  present  formulation,  however,  there  is  no  need 
to  change  the  mean  temperature  of  the  atmosphere  and  the 
interpretation  just  discussed  is  justified. 
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Figure  19  shows  the  result  of  subtracting  the  sen- 
sible heat  term  from  prescribed  radiation  and  latent  heat 
contributions  in  the  ocean  model.   The  curve  has  the 
general  appearance  of  the  S  '  curve  in  Figure  14,   But  a 
"hump"  appears  in  the  same  area  as  the  "dip"  in  the  atmos- 
phere curve  and  should  be  expected  for  the  same  reasons. 

The  oceanic  heating  formulation  that  was  arrived 
at  empirically  and  through  experimentation  was  utilized  in 
the  remainder  of  the  numerical  study. 
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Figure  19.   Average  of  heating  in  the  ocean  model 
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V.   PROCEDURES 

A.   PHASE  I 

1 .   Variations  of  the  Combined  Model 

It  became  evident  during  the  short  term  experiments 
discussed  in  the  last  section  that  the  computing  time  re- 
quired to  achieve  a  state  of  quasi-equilibr ium  was  excessive 
Bryan  (1969)  has  shown  through  scale  analysis  that  the  ver- 
tical advective  time  scale  is  on  the  order  of  decades  while 
the  vertical  diffusive  time  scale  is  on  the  order  of  cen- 
turies.  Both  of  these  time  scales  must  be  relevant  in  the 
integration  of  the  coupled  ocean  and  atmosphere  over  long 
periods  of  time.   The  relatively  shallow  ocean  depth  and 
stratified  temperature  structure  specified  as  the  initial 
conditions  enhanced  the  approach  to  an  equilibrium  state. 

The  combined  model  with  a  200  km  grid  distance  and 
time  increments  as  previously  described  required  approxi- 
mately one  hour  of  computing  time  for  every  100  days  of 
simulation  on  the  IBM  360/67  computer.   This  would  require 
365  hours  of  operating  time  for  100  years  of  simulation  in 
this  phase  of  the  experiment,  and  this  was  prohibitive. 

Several  ways  to  reduce  computing  time  requirements 
were  considered.   These  considerations  involved  the  ocean 
model  since  the  atmospheric  model  was  formulated  simply 
and  could  be  integrated  rapidly.   In  a  combined  model  such 
as  this,  almost  any  change  in  either  time  increment  or  grid 
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distance  in  one  part  of  the  model  would  require  an  adjust- 
ment of  that  quantity  in  the  other  portion  of  the  model. 

A  time  reduction  could  be  realized  by  linearizing 
the  horizontal  equations  of  motion,  as  was  done  by  Haney  in 
the  first  phase  of  his  ocean  experiments.   This  change  was 
incorporated  into  the  model,  along  with  the  assumption  that 
the  vertical  mean  motion  was  steady  so  there  was  no  require- 
ment to  relax  for  \p     ,  the  time  rate  of  change  of  the  stream 
function.   It  was  this  feature  that  produced  the  time  reduc- 
tion.  This  simplified  form  of  the  combined  model  was 
tested  in  a  several  hundred  day  simulation  and  the  model 
advanced  140  days  after  one  hour  of  computing  time.   While 
this  was  a  fair  improvement  in  computing  requirements  it 
was  not  enough  to  warrant  changing  the  model  in  this  way. 

Changing  the  distance  between  grid  points  was  con- 
sidered.  The  model  previously  described  had  a  200  km  grid 
distance  in  both  x  and  y  requiring  horizontal  fields  of  775 
points  (25  by  31).   If  the  grid  distance  was  changed  to 
400  km  in  the  y  direction  and  remained  200  km  in  the  x 
direction,  the  computing  time  would  be  reduced.   The  reduc- 
tion was  due  to  the  fewer  number  of  calculations  to  be  made 
(now  400  grid  points)  and  the  doubling  of  the  time  incre- 
ment in  the  atmospheric  portion  of  the  model.   For  predic- 
tion purposes  the  atmosphere  had  grid  points  only  in  the  y 
direction  and  if  this  increment  was  doubled,  the  time  step 
could  also  be  doubled  and  the  solution  still  remain  stable 
in  a  linear  sense.   The  time  step  in  the  ocean,  however, 
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must  remain  the  same  since  the  grid  distance  in  x  was  not 
changed.   This  version  of  the  model  was  tested  and  produced 
a  320  day  forecast  in  one  hour.   This  was  a  substantial 
time  savings  but  it  was  believed  that  an  even  more  acceler- 
ated pace  could  be  obtained. 

The  final  modification  to  the  original  combined 
model  was  made  by  changing  the  grid  distance  in  both  the  x 
and  y  directions  to  400  km.   When  the  time  steps  were 
doubled  to  one  hour  in  the  atmosphere  and  18  hours  in  the 
ocean,  computing  time  was  reduced  by  a  factor  of  about 
seven.   A  test  showed  that  this  version  of  the  model  pro- 
duced a  666  day  simulation  in  one  hour  of  computing  time. 
This  was  a  satisfactory   improvement  and  this  form  of  the 
model  was  utilized  in  the  early  stages  of  the  experiment. 
2.   The  Coarse  Grid  (400  km)  Model 

Munk  (1950)  has  shown  that  the  development  of  a 
western  (frictional)  boundary  current  of  several  hundred 
kilometers  is  to  be  expected.   It  was  assumed  that  the 
western  boundary  current  was  the  width  of  the  new  grid 
distance,  400  km.   It  has  been  shown  by  Munk  that  the  width 
of  the  boundary  current  (L  )  is  proportional  to  the  cube 
root  of  the  ratio  of  the  lateral  momentum  diffusion  coeffi- 
cient to  the  gradient  of  Coriolis  parameter.   This  can  be 
expressed  as 

\  C  Ab1/3    •  cs-i) 
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The  ratio  of  the  western  boundary  width  to  the  grid  size, 

3,1/3 

L  /Ax,  is  proportional  to  IA  /3(Ax)~]    .   In  order  to 

keep  this  ratio  the  same  in  both  the  400  km  and  200  km  grid 
models,  and  thereby  have  similar  horizontal  resolution  in 
both  models,  the  following  must  hold, 


M 


=  8  A 


(400) 


M 


(200) 


(5.2) 


o     o       _  -I 

Since  a  value  of  1  x  10   cm   sec    was  used  for  A,  with  the 


200  km  grid,  A 


M 


x  10   was  used  as  an  initial  value  in 


this  model.   Similarly  A   was  increased  so  that  a  value  of 

8    2     -1 
A   =  1  x  10   cm   sec    was  used  in  the  coarse  grid  model. 
H 

The  initial  conditions  discussed  previously  were 
modified  to  conform  to  the  new  13  by  16  point  grid  and  the 
long  term  integration  was  started.   Because  of  the  large 
horizontal  eddy  diffusion  of  heat,  the  meridional  sea  sur- 
face temperature  gradient  dropped  from  its  initial  value  of 
22  degrees  over  the  6000  km  length  to  18  degrees  in  the 
first  400  days.   To  have  an  active  air-sea  interaction 
model  the  meridional  temperature  gradient  in  the  ocean  must 
be  large  enough  to  force  strong  meridional  temperature 
gradients  in  the  atmosphere  so  that  a  thermal  wind  larger 
than  that  required  for  baroclinic  instability  can  exist. 
Such  a  large  temperature  gradient  could  be  maintained  in 
the  ocean  only  by  reducing  the  diffusion  coefficient.   How- 
ever, that  coefficient  cannot  be  arbitrarily  reduced  because 
in  the  reduction,  the  western  boundary  width  as  given  by 
(5.1)  may  become  smaller  than  the  grid  size. 
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A  stepwise  reduction  in  the  lateral  diffusion  co- 
efficient was  made.   The  solution  was  examined  for  computa- 
tional stability  after  each  step.   At  the  end  of  the  first 
400  days,  the  lateral  diffusion  coefficients  were  halved, 
the  decrease  in  meridional  temperature  gradient  ceased  and 
the  solution  remained  stable.   At  the  end  of  1200  days  the 
coefficients  were  again  reduced  and  at  the  end  of  7200  days 

o 

they  were  reduced  a  final  time  to  values  of  A.,  =  1  x  10 

M 

o 

and  A   =  .25  x  10   with  no  effect  on  the  computational 
H 

stability.   It  was  evident,  however,  that  with  the  400  km 
grid,  the  western  boundary  current  was  only  marginally  re- 
solved . 

The  model  equations  were  integrated  over  the  coarse 
grid  for  a  period  of  20,200  days  (~55  years).   The  volume 
average  of  temperature  in  the  ocean  basin  was  one  parameter 
used  to  measure  the  approach  to  equilibrium.   This  water 
temperature  had  slowly  cooled  from  an  initial  value  of 
6.3076  degrees  to  5.2017  degrees  at  day  20,200.   In  the 
early  stage  of  adjustment  the  cooling  was  at  a  very  rapid 
rate,  5.2  degrees  per  century,  but  by  the  end  of  this  part 
of  the  experiment,  the  net  cooling  rate  had  reached  zero. 

Another  parameter  monitored  was  the  meridional 
gradient  of  sea  surface  temperature.   This  quantity  had 
stabilized  at  21.6  degrees  per  6000  km  after  day  15,000 
and  remained  so  until  the  end  of  this  time  period.   During 
this  time  the  atmosphere  remained  very  active.   A  thermal 
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wind, (U-  -  U  /2),  with  maximum  values  greater  than  11  m 
sec   ,  which  is  well  above  the  critical  value  for  baroclinic 
instability,  was  present  in  the  mid-portion  of  the  domain. 
A  strong  upper  level  jet  and  weaker  lower  level  jet  were 
present  at  the  same  latitude.   The  state  of  the  atmosphere 
and  ocean  was  very  encouraging  and  the  data  for  day  20,200 
was  used  to  initialize  the  fine  grid  (200  km)  model  to  con- 
tinue the  experiment.   More  discussion  of  the  state  of  the 
atmosphere-ocean  system  at  that  time  is  included  in  the 
section  on  discussion  of  results. 
3 .   Data  Retention 

Prior  to  continuing  with  the  description  of  the  pro- 
cedure used  with  the  fine  grid  model,  a  brief  discussion 
of  data  retention  must  be  included.   In  a  project  requiring 
large  amounts  of  computing  time,  it  is  mandatory  that  a 
great  deal  of  consideration  be  given  to  which  data  should 
be  saved.   One  must  decide  not  only  what  parameters  to 
retain,  but  how  often  and  in  what  form.   It  is  virtually 
impossible  to  repeat  an  experiment  to  extract  parameters 
that  were  not  saved. 

The  data  obviously  must  be  stored  for  analysis  at 
the  end  of  the  experiment.   Data  must  also  be  retained  for 
use  in  restarting  the  model.   It  should  not  be  expected 
that  a  long-term  integration  of  this  nature  could  be  com- 
pleted in  one  computing  time  period.   The  equations  must 
be  integrated  and  data  stored  periodically  so  that  after  a 
reasonable  operating  time,  such  as  one  hour,  the  model  can 
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be  stopped  and  restarted  at  a  later  time  using  the  last 
data  in  storage  as  initial  values. 

Throughout  the  experiment  with  the  coarse  grid 
model,  data  were  stored  at  100  day  intervals.   The  ocean 
variables  saved  were  the  stream  function  from  which  the 
vertical  mean  current  was  obtained,  the  zonal  vertical 
shear  current,  u',  the  meridional  vertical  shear  current, 
v',  and  the  temperature,  T,  at  all  three  levels.   The  at- 
mospheric variables  saved  were  the  Fourier  coefficients, 
A  through  F.   The  retention  of  such  a  large  quantity  of 
data  posed  a  problem,  but  it  was  remedied  through  the  judi- 
cious use  of  magnetic  tapes.   As  one  section  of  the  experi- 
ment was  completed,  the  data  that  would  not  be  needed  until 
the  final  analysis  was  moved  from  disk  pack  temporary 
storage  to  a  set  of  magnetic  tapes  where  it  was  permanently 
stored.   The  necessity  for  this  became  even  more  evident  in 
the  fine  grid  model  integrations  when  the  fields  were  in- 
creased from  208  data  points  per  level  to  775  points  per 
level.   Saving  this  much  data  would  have  been  impossible 
without  a  large  storage  capability. 

The  last  problem  concerning  data  storage  arose  dur- 
ing the  final  part  of  the  long  term  integration.   After  a 
state  of  quasi-equilibrium  was  approached,  five  points  were 
chosen  for  retention  of  daily  data  for  the  last  six  years 
of  the  experiment  and  were  carefully  examined  by  spectral 
techniques.   Data  from  an  ocean-only  and  an  atmosphere-only 
experiment  were  also  saved  over  the  last  six  years. 
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4 .   Fine  Grid  (200  km)  Model 

The  data  for  day  20,200  from  the  coarse  grid  model 
was  linearly  interpolated  in  both  x  and  y  to  obtain  initial 
fields  for  the  fine  grid  model.   The  lateral  diffusion  co- 
efficients for  both  momentum  and  heat  were  held  at  the 
final  values  as  used  in  the  coarse  grid  model.   While  these 
coefficients  might  have  been  reduced  over  a  time  period  as 
they  were  in  the  coarse  grid,  this  was  not  done  in  order  to 
avoid  the  readjustment  that  would  be  required.   A  further 
change  in  the  diffusion  coefficient  would  increase  the  ad- 
justment period.   It  has  been  shown  by  Manabe  (1969a), 
Gates  (1968)  and  others  that  the  values  chosen  previously 
are  appropriate  for  the  grid  size  used. 

The  equations  were  integrated  over  the  fine  grid 
model  for  another  6500  days  with  data  retained  every  100 
days.   At  that  time,  73  years,  the  ocean-atmosphere  model 
was  approaching  an  equilibrium  state.   The  volume  average 
of  the  temperature  initially  started  to  slowly  decrease 
but  after  a  few  hundred  days  it  commenced  an  increasing 
trend  that  it  maintained  throughout  the  remainder  of  the 
experiment.   The  warming  trend  began  at  the  rate  of  0.8 
degrees  per  century  but  this  decreased  to  about  0.5  degrees 
per  century  at  73  years.   Compared  to  the  warming  rate  of 
2  degrees  per  century  arrived  at  by  the  ocean  in  Manabe 
and  Bryan's  combined  model,  this  was  considered  an  accept- 
able time  to  start  the  last  six  years  of  the  combined 
model  experiment. 


96 


Beginning  at  the  73rd  year,  a  final  2200  day  (6 
years)  simulation  was  produced  with  the  fine  grid  model. 
It  was  during  this  time  that  not  only  were  all  the  data 
recorded  every  100  days  but  for  five  selected  points,  the 
data  were  recorded  every  1.125  days  (every  three  oceanic 
time  steps).   These  points  were  selected  as  being  represen- 
tative of  features  that  existed  in  various  parts  of  the 
ocean.   Figure  20  shows  the  location  of  the  five  points, 


72°N 


Figure  20.   Location  of  points  where  selected  data  were 
recorded  every  1.125  days. 
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lettered  A  through  E.   Point  A  is  the  northernmost  point  in 
the  eastern  sector  where  sea  surface  temperatures  were 
warmer  than  the  initial  values  and  where  strong  sinking 
motion  had  occurred.   Point  B  is  in  the  center  of  the  domain 
where  the  cyclonic  scale  atmospheric  disturbances  had  their 
greatest  effect.   Point  C  is  in  the  western  boundary  current 
and,  if  the  model  domain  is  likened  to  the  North  Atlantic 
Ocean,  this  point  is  close  enough  to  an  observation  site 
used  by  Woods  Hole  Oceanographic  Institution  so  that  some 
comparisons  can  be  made.   Point  D  is  an  area  in  the  south 
central  portion  of  the  domain  and  would  be  equivalent  to 
an  area  near  Bermuda.   Point  E  is  in  the  warm  core  of  the 
ocean  close  to  the  southwestern  corner  of  the  domain.   The 
data  at  these  points  were  examined  spectrally  and  the  dis- 
cussion of  this  technique  is  included  in  a  later  section. 
At  the  conclusion  of  this  part  of  the  experiment, 
the  combined  model  had  been  used  for  a  79  year  integration 
of  the  equations.   This  completed  phase  I  of  this  project. 

B.   PHASE  II 

This  phase  of  the  experiment  consisted  of  integrating 
the  equations  of  the  ocean-only  model  for  a  six-year  period 
using  the  fine  grid.   To  simulate  the  coincidence  of  this 
time  period  with  the  final  six  years  of  the  combined  model, 
data  from  the  combined  model  were  used  for  initialization. 
Initial  fields  of  temperature  and  velocities  were  taken 
from  day  26,700. 
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A  fixed  surface  stress  and  an  atmospheric  temperature 
had  to  be  imposed.   Rather  than  use  climatology,  the  input 
stress  and  temperature  were  obtained  from  the  combined 
model.   In  preparation  for  this  phase  of  the  experiment, 
the  surface  stress,  and  atmospheric  surface  air  temperature 
were  averaged  over  a  100  day  period  centered  at  day  26,700. 
These  average  values  were  then  held  fixed  during  the  entire 
time  integration  of  the  ocean-only  model. 

Data  were  retained  in  exactly  the  same  manner  as  in  the 
combined  model.   Full  fields  of  the  ocean  variables  were 
stored  every  100  days  and  data  for  the  five  special  points 
were  stored  every  1.125  days.   The  variables  were  examined 
and  comparisons  were  made  with  the  results  of  Phase  I. 
These  results  are  discussed  in  a  later  part  of  this  work. 

C.   PHASE  III 

The  equations  for  the  atmosphere- only  model  were  inte- 
grated over  the  same  six  year  time  span  with  the  fine  grid. 
Initial  values  of  A  through  F  were  taken  from  day  26,700 
of  the  combined  model  data  and  the  100  day  average  sea  sur- 
face temperatures  were  imposed  as  a  lower  boundary  condi- 
tion.  Data  were  saved  at  the  same  times  as  in  the  first 
two  phases . 

Calculations  of  energy  and  energy  transformations  were 
made  and  are  described  later.   The  data  from  the  five  points 
were  analyzed  spectrally  and  this  too  is  presented  ii\  a 
later  section.   With  the  completion  of  this  phase  of  the 
experiment,  the  long-term  integrations  were  concluded. 
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D.   REVIEW  OF  COMPUTING  TIME  REQUIREMENTS 

Since  the  integrations  were  conducted  with  a  variety  of 
numerical  models,  it  is  germain  at  this  point  to  review  the 
computing  time  requirements  for  each  model.  This  review  is 
presented  in  tabular  form  in  Table  2. 


MODEL 

DESCRIPTION 

LENGTH  OF  SIMULATION 
FOR  1  HR  COMPUTING  TIME 

Combined 
(Fine  Grid) 

Ax=Ay=200km 

At    =30  min 

atm 

At    =9  hours 
ocn 

101  days 

Combined 
(Coarse  Grid) 

Ax=Ay=400km 

At    =1  hour 
atm 

At    =18  hours 
ocn 

666  days 

Ocean  Only 
(Fine  Grid) 

Ax=Ay=200km 
At=  9  hours 

308  days 

Atmosphere  Only 
(Fine  Grid) 

Ay=  200km 
At=   30  min 

697  days 

Combined 
(Linear) 
(Not  used) 

Ax=Ay=200km 

At    =  30  min 
atm 

At    =9  hours 
ocn 

linear  momentum 

equations 

140  days 

Combined 
(Mixed  Grid) 
(Not  used) 

Ax=200km 

Ay=400km 

At    =1  hour 
atm 

At    =9  hours 
ocn 

320  days 

TABLE  2.   Review  of  computing  time  requirements  for 
variations  of  the  numerical  models  used  in  this  experiment 
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It  can  be  readily  deduced  from  these  time  figures  that 
the  coupling  processes  required  lengthy  calculations.   The 
method  by  which  the  sea  surface  temperatures  are  transferred 
into  wave  number  space  required  many  time  consuming  opera- 
tions for  every  oceanic  time  step.   The  other  summations 
and  averaging  procedures  all  contributed  to  a  joining  pro- 
cess that  caused  the  combined  model  to  run  significantly 
slower  than  either  of  its  individual  components. 
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VI.   DISCUSSION 

The  data  that  was  collected  during  the  long  term  inte- 
gration can  be  placed  in  two  categories.   The  full  fields 
of  data  and  the  data  collected  at  the  five  designated 
points  differ  widely  in  frequency  of  observations  and  area 
of  coverage.   Because  of  this,  the  methods  by  which  the 
data  were  to  be  examined  differ  considerably.   The  full 
fields  of  data  that  were  retained  at  100  day  intervals  were 
examined  in  a  synop tic- type  study,  while  the  five-point 
data,  recorded  at  daily  intervals,  were  examined  using 
spectral  techniques. 

In  the  synoptic  study  the  analyzed  data  were  compared 
not  only  to  observed  climatological  data,  but  also  to  the 
results  obtained  from  the  various  phases  of  this  numerical 
experiment.   The  atmospheric  energy  balance  was  examined 
and  comparisons  were  made  between  phase  I  and  phase  III  and 
with  the  work  discussed  in  Chapter  II.   In  the  spectral 
study, power  spectrum  estimates  that  were  obtained  from  data 
at  one  particular  point  were  compared  with  estimates  ob- 
tained from  data  at  other  points  during  the  same  phase  of 
the  experiment.   Similarly,  spectral  estimates  of  one 
variable  were  compared  with  estimates  of  the  same  variable 
at  the  same  point  in  a  different  phase  of  the  work.   Some 
comparisons  were  made  with  spectral  estimates  calculated 
from  observed  data  collected  by  the  Woods  Hole  Oceanographic 
Institution . 
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A.   SYNOPTIC  STUDIES 

Throughout  the  experiment  full  fields  of  the  stream  func- 
tion, vertical  shear  current  components,  and  temperatures 
were  retained  as  oceanic  data  and  the  Fourier  coefficients 
were  saved  as  the  atmospheric  data.   These  variables,  along 
with  fields  of  oceanic  vertical  motion  and  atmospheric  zonal 
velocity,  vertical  motion,  mean  meridional  velocity  and 
disturbance  amplitude  and  phase  information  were  printed 
every  50  or  100  days.   Many  of  these  variables  were  examined 
in  detail  and  discussion  of  the  results  is  contained  in  the 
subsequent  sections. 

1 .   The  400  km  Grid  Combined  Model  Experiment 

The  400  km  grid  combined  model  was  merely  a  tool  that 
was  used  to  expeditiously  integrate  the  model  equations  to 
reach  a  state  of  quasi-equilibr ium .   The  ocean-atmosphere 
interaction  at  the  end  of  this  part  of  phase  I  was  not  of 
primary  interest.   The  state  of  the  system  at  the  end  of 
day  20,200  (55  years)  was  to  be  used  only  to  initialize  the 
fine  grid  model  from  which  more  realistic  and  more  meaning- 
ful results  may  be  obtained.   However,  a  brief  discussion 
of  the  conditions  to  which  this  form  of  the  combined  model 
evolved  is  necessary  to  explain  why  the  2  00  km  grid  was 
needed  to  complete  this  phase. 

The  initial  conditions  for  this  part  of  the  experi- 
ment have  previously  been  described.   The  initial  adjustment 
of  sea  surface  temperatures  was  due  to  the  large  horizontal 
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eddy  diffusion  coefficient  then  in  use.   The  reduction  of 
this  coefficient  and  the  necessity  for  maintenance  of  a 
large  meridional  sea-surface  temperature  gradient  was  al- 
ready discussed.   The  coefficient  of  lateral  diffusion 
could  not  be  made  small  enough  with  the  400  km  grid  so  that 
nonlinear  terms  would  be  important.   Bryan  and  Cox  (1967) 
have  defined  an  "effective"  Reynolds  number,  Re,  as  the 
Reynolds  number  associated  with  the  lateral  turbulent  mix- 
ing.  The  nonlinearity  of  the  solution  is  dependent  on  the 
size  of  the  Reynolds  number  defined  as  follows: 


Re  E 


U*L 


(6.1) 


where  U*  is  the  velocity  component  depending  on  the  density 

gradient  , 

L  is  a  scale  length,  and 

A,,  is  the  coefficient  of  lateral  momentum  diffusion. 
M 

Bryan  and  Cox  conducted  a  series  of  experiments 
with  values  of  Re  =  0.8,  which  gave  a  highly  diffusive, 
linear  solution  and  Re  =  8.0,  which  yielded  a  moderately 
diffuse,  slightly  nonlinear  solution.   The  maximum  value 
of  Re  used  in  this  experiment  was  6.5  and  it  is  believed 
that  a  value  greater  than  10  was  necessary  to  achieve  a 
truly  nonlinear  solution.   It  is  easily  shown  that  with  the 

400  km  grid,  continued  reduction  of  A„  produces  a  situation 

M 

where  the  linear  western  boundary  current  becomes  unresolv- 
able  [see  Equation  (4.1)]  prior  to  the  occurrence  of  non- 
linear effects.   Therefore,  the  solution  with  the  400  km 
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grid  is  similar  to  the  viscous  linear  solution  described 
by  Munk  (1950).   Figure  21  shows  the  ocean  stream  function 
obtained  at  day  20,200  (about  55  years).   According  to 
Munk,  the  streamlines  should  represent  an  exponentially 
decaying  series  of  vorticies  centered  on  the  latitudinal 
axis.   The  vorticies  are  present  across  the  width  of  the 
domain  but  do  not  appear  to  decay  as  expected.   Truncation 
error  is  a  probable  cause  for  insufficient  decay  in  this 
area.   In  the  gyre  north  of  45N,  the  maximum  value  of  the 
vertically  averaged  meridional  current  was  on  the  order  of 
10  cm  sec   .   In  the  southern  gyre  the  maximum  value  of  the 
same  quantity  was  about  15  cm  sec   .   The  magnitude  of  the 
total  (vertical  shear  plus  vertical  mean  parts)  meridional 
currents  in  the  northern  and  southern  gyres  were  10  and  35 
cm  sec    respectively.   The  higher  latitudes  in  the  model 
were  more  barotropic  and  the  vertical  mean  current  domina- 
ted, while  the  vertical  shear  current  dominated  in  the 
highly  baroclinic  areas  of  the  central  portion.   In  the 
north  the  vertical  shear  current  was  about  zero,  while  in 
the  south  the  vertical  shear  current  was  about  20  cm  sec 
directed  northward. 

Figure  22  shows  the  sea  surface  temperature  field 
at  day  20,200.   The  mean  meridional  temperature  gradient 
over  the  latitudinal  domain  was  21.4  degrees  per  6000  km. 
The  temperatures  at  both  the  north  and  south  boundaries  had 
decreased  from  their  initial  values.   Upwelling  along  the 
central  portion  of  the  northern  boundary  was  responsible 
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Long  i  f  ude 


Figure  21.   Ocean  stream  function  obtained  at  day  20,200 
from  400  km  grid  combined  model. 
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L  ong  i  f  u  d  e 


Figure  22.   Sea  surface  temperature  at  day  20,200,  (C). 
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for  this  cooling  at  the  surface.   The  small  scale  spatial 
fluctuation  in  the  temperature  field  in  low  latitudes  ap- 
pears to  be  a  result  of  the  linear  viscous  solution  in 
which  there  is  alternating  northerly  and  southerly  flow  and 
therefore  northerly  and  southerly  advection  of  heat. 

Figure  23  shows  the  zonal  winds  at  three  levels  in 
the  atmosphere.   The  curve  labeled  (a)  depicts  the  zonal 
mean  surface  wind.   The  easterly,  westerly,  easterly  pattern 
from  north  to  south  is  qualitatively  like  the  patterns  in 
the  real  atmosphere.   However,  the  latitudes  at  which  the 
direction  of  flow  change  occur  are  slightly  north  of  the 
actual  latitude  of  change.   Curve  (b)  is  the  zonal  wind  at 
750mb  and  curve  (c)  is  that  quantity  at  the  250  mb  level. 
The  increase  of  westerlies  with  height  in  the  central 
latitudes  is  in  thermal  wind  balance  with  a  temperature 
structure  of  warm  air  in  the  south  and  cold  air  in  the 
north . 

While  the  structure  of  the  atmosphere-ocean  system 
that  evolved  from  integrating  the  model  equations  over  the 
400  km  grid  was  not  extremely  accurate,  it  did  suffice  for 
initializing  the  integration  of  the  fine  grid  model. 
2 .   The  200  km  Grid  Combined  Model  Experiment 

The  fine  grid  or  200  km  grid  part  of  phase  I  was 
initialized  with  data  from  day  20,200  of  the  400  km  grid 
model  after  spatial  interpolation.   The  model  equations 
were  then  integrated  in  time  for  8700  days  (approximately 
24  years).   The  adjustment  to  the  fine  grid  was  evident  in 
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Figure  23.   Zonal  atmospheric  winds  in  m  sec    from  400  km 
grid  model  at  a)  surface;  b)  750  mb ;  c)  250  mb . 
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a  slight  initial  decrease  of  the  oceans  volume  averaged  tem- 
perature.  After  a  few  hundred  days,  this  temperature  start- 
ed increasing  and  at  day  21,000  the  sea  was  warming  at  a 
rate  of  0.8  degrees  per  century.   This  warming  rate  was  main' 
tained  for  over  1000  days  and  then  the  rate  decreased  grad- 
ually to  0.5  degrees  per  century  at  day  28,900,  the  last  day 
of  the  experiment.   This  warming  rate  was  caused  by  a  down- 
ward heat  flux  across  the  air-ocean  interface  of  about  2.3 

-2    -1 
cal  cm   day   .   While  the  experiment  could  have  been  con- 
tinued until  no  heat  flux  was  observed,  it  was  believed  that 
the  large  amount  of  computing  time  required  to  accomplish 
this  would  not  have  been  worth  the  slight  change  in  the 
final  results . 

At  the  end  of  the  400  km  grid  model  experiment,  the 

total  disturbance  energy  in  the  atmosphere  was  statistical- 

2     -2 
ly  steady,  oscillating  about  50  m   sec   .   During  the  first 

few  hundred  days  of  the  fine  grid  integrations  the  average 

2     -2 
disturbance  energy  jumped  to  about  70  m   sec   ,  but  then 

quickly  returned  to  the  statistically  steady  state  which 

had  previously  been  observed.   The  initial  cooling  of  the 

ocean  imparted  an  added  amount  of  heat  to  the  atmosphere 

in  the  southern  half  of  the  region  generating  more  zonal 

potential  energy  which  the  baroclinic  processes  converted 

into  disturbance  energy.   After  the  cooling  ended,  the 

warming  of  the  ocean  apparently  occurred  rather  uniformly 

in  space  resulting  in  the  generation  of  about  the  same 

amount  of  available  potential  energy  as  had  been  generated 
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prior  to  reducing  the  grid  size.   More  discussion  of  the 
energy  balance  in  the  fine  grid  model  is  in  a  later  section. 

The  mean  state  of  the  atmosphere  over  the  last  100 
days  of  this  phase  of  the  experiment  is  shown  in  Figures  2  4 
to  26.   These  maps  along  with  Figures  28  to  30  that  will  be 
discussed  shortly  are  deviation  values  of  the  height  fields 
at  various  levels  obtained  from  the  geostrophic  condition. 
The  zero  contour  and  plus  and  minus  deviation  contours  are 
drawn.   These  maps  are  displayed  in  this  way  because  the 
deviation  of  the  stream  function  is  predicted  by  the  quasi- 
geostrophic  atmospheric  equations.   The  quantity  ip  appears 
as  Vip  in  all  terras,  and  can  be  represented  as  V(ip'(x,y,t) 
+  \\)  ( t )  )  ,  a  deviation  and  mean  part  of  the  stream  function. 
Therefore,  any  mean  ip  could  be  added  to  (or  subtracted 
from)  the  deviation  quantity. 

Figure  24  shows  the  extrapolated  distribution  of 
the  1000  mb  height  deviation  analyzed  at  12-meter  intervals 
There  is  a  semi-permanent  mid-latitude  cyclone  over  the 
eastern  ocean  and  a  high  at  the  same  longitude  south  of 
the  low.   A  high  is  present  in  the  northernmost  section  of 
the  region.   This  is  worthy  of  note  since  there  is  no  cold 
land  mass  modeled  into  the  system  to  simulate  a  "Greenland" 
type  effect.   Figure  25  shows  the  distribution  of  the  250 
mb  height  deviation  analyzed  at  60-meter  intervals.   The 
flow  is  zonal  throughout  the  central  portion  of  the  domain 
with  a  cyclonic  polar  vortex  present  in  the  north  and  a 
weak  ridge  present  at  the  same  longitude  in  the  south. 
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Figure  24.   Distribution  of  100  day  mean  1000-mb  height 
deviation  in  phase  I . 
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Figure  25.   Distribution  of  100  day  mean  250  mb  height 
deviation  in  phase  I . 
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Figure  26  shows  the  distribution  of  the  750  mb  height  de- 
viation analyzed  at  30-meter  intervals.   The  flow  is  again 
zonal  throughout  the  central  portion  of  the  region.   How- 
ever, the  lowest  height  is  now  in  the  northeastern  corner 
of  the  area  and  the  low  latitude  high  is  present  in  only 
the  southernmost  part  of  the  domain. 

The  vertical  structure  of  the  standing  disturbance 
is  such  that  the  1000  mb  high  in  the  polar  regions  lies 
below  the  250  mb  low  the  same  as  in  the  real  atmosphere. 
The  warm  low  pressure  area  at  the  southern  boundary  at  the 
1000  mb  level  becomes  an  area  of  high  pressure  at  the  250 
mb  level. 

The  zonal  wind  averaged  over  the  same  100  day  per- 
iod is  shown  in  Figure  27.   The  line  labeled  (c)  shows  a 
westerly  flow  at  all  latitudes  at  the  250  mb  level  with  the 
jet  located  at  45N,  coinciding  with  the  strong  height 
gradient  in  that  region  as  seen  in  Figure  24.   Curve  (b) 
shows  the  weak  polar  easterlies  at  750  mb  level  to  the 
north  of  the  low  at  60N.   The  remainder  of  the  flow  is 
westerly  with  a  weak  jet  at  about  45N.   Line  (a)  shows  the 
zonal  flow  at  the  1000  mb  level.   The  polar  easterlies 
extend  from  the  northern  boundary  to  a  latitude  just  north 
of  the  cyclonic  activity.   Westerlies  extend  from  that 
point  to  the  latitude  of  the  1000  mb  highs  and  easterlies 
exist  to  the  southern  boundary.   A  comparison  of  Figure  27 
with  Figure  23  shows  that  the  mean  zonal  winds  did  not 
change  significantly  after  reducing  the  grid  size.   This 
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Figure  26.   Distribution  of  100  day  mean  750  mb  height 
deviation  in  phase  I. 
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indicates  that  for  the  atmospheric  model,  a  400  km  grid 
would  have  sufficed  with  little  loss  of  resolution. 

Figure  28  shows  the  distribution  of  the  1000  mb 
height  contours  at  25  meter  intervals  at  day  28,900.   A  lo 
dominates  the  western  half  of  the  region  at  this  time  and  a 
high  covers  the  eastern  half.   This  figure  was  compared  to 
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Figure  27.   Zonal  atmospheric  winds  in  m  sec 
grid  model  at  a)  1000  mb ;  b)  750  mb ;  c)  250  mb . 


the  series  of  1000  mb  figures  presented  by  Phillips  (1956) 
and  the  gross  features  are  almost  identical  to  the  early 
stages  of  his  model  atmospheric  development.   The  main  sur- 
face troughs  and  ridges  are  orientated  so  as  to  lean  back 
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toward  the  northwest  in  the  northern  half  and  slightly  to- 
ward the  southwest  in  the  southern  half  of  the  region.  The 
convergence  of  momentum  in  the  central  portion  increases 
the  westerlies  and  accounts  for  the  mid-latitude  jet  ob- 
served in  this  area  (Figure  29). 

Figure  29  shows  the  distribution  of  the  250  mb 
height  contours  at  60  meter  intervals  for  day  28,900  (79 
years).   A  wave  of  moderate  amplitude  is  superimposed  on 
the  zonal  flow  with  a  trough  present  over  the  western  half 
of  the  region  and  a  ridge  dominating  the  eastern  half.   The 
polar  vortex  observed  in  Figure  25  exists  at  this  time 
also.   The  single  wave  structure  of  the  atmospheric  model 
forces  a  pattern  similar  to  this  to  be  found  at  any  given 
instant.   The  disturbance  moves  and  its  amplitude  varies 
with  time  from  a  zonal  pattern  to  a  highly  amplified  wave 

pattern.   At  this  time  the  disturbance  kinetic  energy  had  a 

2     -2 
value  of  54  m   sec   .   This  is  quite  large  compared  to  the 

disturbance  kinetic  energy  of  the  more  zonal  flow  in  Figure 

25.   However,  at  other  times  in  the  experiment  higher 

values  of  this  quantity  were  observed.   It  is  believed  that 

the  atmospheric  state  shown  in  Figure  28  is  representative 

of  a  flow  in  which  a  wave  disturbance  of  average  amplitude 

is  imbedded.   At  times  the  amplitude  of  the  wave  will  be 

greater  than  it  is  at  day  28,900  and  at  other  times  less. 

Also,  the  wave  will  propagate  across  the  domain  and  yield 

a  time  mean  upper  level  flow  as  shown  in  Figure  25. 
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Figure  28.   Distribution  of  1000  mb  height  deviation  at  day 
28,900  in  phase  I. 
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Figure  29.   Distribution  of  250  mb  height  deviation  at  day 
28,900  in  phase  I. 
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Figure  30,  the  distribution  of  the  750  mb  height 
contour  at  30-meter  intervals  at  day  28,900,  shows  a  closed 
cyclonic  circulation  directly  beneath  the  trough  at  250  mb . 
The  amplitude  of  the  wave  motion  is  about  the  same  as  at 
250  mb  and  a  closed  high  is  present  beneath  the  upper  level 
ridge  . 

The  suggestion  that  the  single  wave  in  the  x-direc- 
tion  would  produce  a  series  of  cyclonic  scale  disturbances 
moving  over  the  domain  of  the  atmosphere-ocean  model  was 
born  out  by  these  figures  plus  examination  of  preceding 
maps.   The  wave  moves  eastward  with  a  speed  of  about  600  km 
per  day.   This  is  considerably  slower  than  the  wave  motion 
in  Phillips'  experiments. 

The  vertical  heat  flux  was  discussed  earlier  as  be- 
ing proportional  to  the  atmosphere-ocean  temperature  differ- 
ence.  However,  due  to  the  method  by  which  the  sensible 
heat  exchange  is  calculated,  the  vertical  heat  flux  is 
actually  proportional  to  the  difference  between  the  atmos- 
pheric  surface  temperature  and  the  Fourier  analyzed  sea 
surface  temperature.   This  temperature  difference,  averaged 
over  the  last  100  days  of  the  experiment,  is  shown  in 
Figure  31.   In  general,  there  is  cooling  of  the  atmosphere 
in  the  western  sector  and  warming  of  the  atmosphere  in  the 
east.   While  this  heating  may  be  representative  of  what 
occurs  in  the  northernmost  latitudes  of  the  real  atmosphere, 
it  does  not  accurately  portray  what  occurs  throughout  the 
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Figure  30.   Distribution  of  750  mb  height  deviation  at  day 
28,900  in  phase  I. 
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Figure  31.   100  day  average  difference  between  the  Fourier 
transformed  sea  surface  temperature  and  the  atmosphere  sur> 
face  temperature  computed  in  phase  I.   Heating  of  the 
atmosphere  occurs  where  values  are  positive. 
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majority  of  the  oceanic  regions  where  cold  atmospheric  air 
masses  move  from  over  a  cold  continent  to  a  position  over 
warmer  water  and  are  heated  from  below.   The  lack  of  colder 
continents  in  the  model  and  the  use  of  cyclic  boundary  con- 
ditions makes  it  impossible  to  simulate  this  phenomena  in 
this  simple  model.   The  use  of  the  Fourier  analyzed  sea 
surface  temperature  also  contributes  to  the  unrealistic 
vertical  heat  flux  pattern  of  the  model.   The  strong  tempera- 
ture gradient  observed  in  the  real  ocean  along  the  east 
coast  of  continents  and  observed  in  the  model  along  the 
western  boundary  of  the  ocean  in  Figure  34,  had  little 
effect  on  the  heat  flux  pattern  that  evolved.   The  Fourier 
analysis  smoothed  this  gradient  and  although  the  gradient 
did  influence  the  new  sea  surface  temperature  somewhat,  it 
did  not  have  the  same  influence  observed  in  the  real  atmos- 
phere-ocean system.   The  Fourier  analyzed  sea  surface  tem- 
perature was  used  in  the  heat  exchange  to  be  consistent 
with  the  simple  one-wave  formulation  of  the  atmosphere. 
The  cyclic  boundary  condition  forces  the  air  to  pass  from 
a  region  of  warm  water  on  the  eastern  boundary  immediately 
to  the  region  of  cooler  water  in  the  western  area.   The 
overall  effect  that  the  resulting  heat  flux  pattern  had  on 
the  results  of  this  experiment  are  hard  to  determine.   The 
strong  heating  observed  in  nature  on  the  western  side  of  the 
ocean  in  the  mid-latitudes  does  not  occur.   The  overall 
effect  may  be  to  smooth  out  the  effect  of  the  heating  on 
the  atmospheric  temperature  disturbance  (only  disturbance 
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quantities  are  of  interest  since  domain  mean  temperatures 
of  the  atmosphere  do  not  change)  and  cause  a  more  zonal 
pattern.   Figure  32  shows  the  time  mean  atmospheric  surface 
temperature  averaged  over  the  last  100  days  of  the  combined 
model  experiment.   The  temperature  in  this  figure  has  a 
more  zonal  appearance  than  was  expected  after  examining  the 
mean  1000  mb  height  deviations  (Figure  24) . 

Figure  33  is  the  100  day  time  mean  of  the  stream 
function  governing  the  vertical  mean  currents.   The  sub- 
arctic gyre,  visible  at  about  55N  in  this  figure,  produces 
southward  flow  along  the  western  boundary  at  that  latitude 
and  northward  return  flow  4-600  km  east  of  the  boundary. 
The  subtropical  circulation  at  about  35N  has  northward  flow 
at  the  western  boundary  and  southward  flow  4-600  km  to  the 
east.   The  convergent  flow  in  the  boundary  causes  a  sharp 
meridional  temperature  gradient  at  about  43N  and  contri- 
butes to  the  low  temperature  near  48N  in  the  western  bound- 
ary . 

Figures  34,  35,  and  36  show  the  100  day  time  means 
of  the  ocean  temperature  at  the  surface,  500  meters  and 
1300  meters,  respectively.   These  maps  are  vertically  con- 
sistent and  the  low  temperatures  just  described  near  48N  can 
be  seen  at  all  three  levels.   The  thermal  ridge  observed 
in  the  northwestern  sector  of  Figure  34,  and  also  in  Figure 
35,  is  caused  by  the  northwestward  flow  on  the  east  side 
of  the  subarctic  gyre  and  somewhat  by  the  upwelling  along 
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Figure  32.   100  day  mean  atmosphere  surface  temperature  in 
phase  I  (C) . 
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Figure  33.   100  day  mean  stream  function  pattern  from 
phase  X  in  units  of  107  cm2  sec""1. 
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Figure  34.   100  day  mean  sea  surface  temperature  (C) 
pattern  in  phase  I. 
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Figure  35.   Temperature  pattern  at  a  depth  of  500  meters  in 
ocean  in  phase  I  (C) . 
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Figure  36.   Temperature  pattern  at  depth  of  1300  meters  in 
ocean  in  phase  1(C). 
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the  western  boundary  seen  in  Figure  37.  which  contains  the 
time  averaged  vertical  motion  field  at  the  depth  of  300 
meters.   In  the  fine  grid,  as  well  as  in  the  coarse  grid 
discussed  earlier,  the  flow  in  the  western  boundary  regions 
in  high  latitudes  is  composed  almost  entirely  of  the  verti- 
cal mean  current.   Here  the  flow  is  of  this  character  and 
moves  across  the  isotherms,  advecting  warmer  water  toward 
the  northwest.   The  large  area  with  horizontally  isothermal 
water  of  about  5C  located  west  of  the  ridging  in  the  north- 
west quadrant  of  the  ocean  domain  can  be  compared  qualita- 
tively to  a  similar  feature  found  east  of  Newfoundland  in 
the  North  Atlantic.   Sverdrup  et  al  (1942)  presented  sea 
surface  temperature  maps  that  contained  this  feature  in 
both  winter  and  summer. 

In  the  northeastern  quadrant  the  flow  is  weak  and 
closely  parallels  the  isotherms.   Here,  in  the  absence  of 
appreciable  vertical  mean  current,  the  vertical  shear  cur- 
rent is  dominant  and  when  the  flow  reaches  the  northeastern 
boundaries  large  sinking  motions  occur  (see  Figure  37)  pro- 
ducing warm  temperature  at  all  depths  there.   Figures  34, 
35  and  36  all  show  warmer  water  in  the  northeast  corner 
due  to  this  phenomenon.   The  isotherm  pattern  produced  by 
the  model  in  this  area  can  be  likened  to  the  real  sea  sur- 
face temperature  pattern  off  the  west  coast  of  the  United 
States.   The  area  of  warmest  surface  water  in  Figure  34  is 
observed  in  the  southwestern  corner  of  the  domain.   This 
area  is  analogous  to  the  Sargasso  Sea  with  its  warm  tem- 
peratures and  weak  currents. 
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Figure  37.   100  day  mean  vertical  motions  in  cm  day    at 
300  meters  in  phase  I. 
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The  temperatures  along  the  western  boundary  are  fur' 
ther  affected  by  the  upwelling  that  takes  place  due  to  the 
thermal  circulation.   The  central  latitude  is  a  region  of 
weak  eastward  flow  in  the  surface  layer  and  westward  flow 
at  the  lower  layers.   The  moderate  upwelling  brings  colder 
water  to  the  surface,  lowering  the  surface  temperature 
along  the  boundary.   An  area  of  weak  sinking  motion  in  the 
south  central  portion  of  the  domain  in  Figure  37  brings 
about  the  warm  8C  water  at  the  500  meter  level  (Figure  35) . 

From  the  stream  function  field  in  Figure  33,  the 
vertical  mean  current  can  be  computed.   The  subarctic  gyre 
has  a  vertical  mean  current  on  the  order  of  20  cm  sec 
The  subtropical  gyre  has  a  vertical  mean  current  on  the 
order  of  15  cm  sec   .   Table  3  shows  the  total  meridional 
velocity,  an  estimate  of  the  vertical  mean  meridional  cur- 
rent as  calculated  from  the  stream  function  and  the  verti- 
cal shear  meridional  current  calculated  as  a  residual,  in 
the  western  boundary  current.   The  vertical  mean  current 
is  almost  entirely  responsible  for  the  meridional  current 
in  the  boundary  layer  associated  with  the  subarctic  gyre, 
while  the  vertical  shear  current  is  dominant  in  the  sub- 
tropical gyre.   The  maximum  value  of  51  cm  sec    in  the 
boundary  current,  representing  a  Gulf  Stream  type  flow,  is 
less  than  that  found  in  the  real  world  oceans  where  100 
cm  sec    have  been  observed  for  this  type  of  a  boundary 
current. 
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Latitude 

Total  Merid. 

Vertical 

Vertical 

Current 

Mean 

Shear 

71.1 

2.6 

2  .8 

-0.2 

67.5 

12.8 

5.6 

7  .2 

63.9 

9.9 

10.0 

-0.1 

60.3 

-0.4 

-1.7 

1.3 

56.7 

-12.0 

-14  .0 

2.0 

53.1 

-20.5 

-19  .0 

-1.5 

49.5 

-18.7 

-14.0 

-4.7 

45.9 

-5.6 

-6.5 

0.9 

42.3 

11.1 

0.8 

10.3 

40.5 

17.7 

3.8 

13.9 

36.9 

37.5 

12  .5 

25  .0 

33.3 

51.1 

15  .5 

35.6 

29.  7 

41.0 

8.0 

33.0 

26  .1 

20.5 

2  .5 

18.0 

22.5 

6.7 

-5  .1 

11.8 

18.9 

1.9 

-0.7 

2  .6 

TABLE  3.   Total  meridional  current  (cm  sec   )  as  produced 
by  the  combined  model;  vertical  mean  current  as  computed 
from  the  model  stream  function;  vertical  shear  current 
calculated  as  a  residual. 
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Figure  38.   100  day  mean  stress  pattern  produced  by  com- 
bined model . 


The  stress  acting  as  a  driving  force  on  the  ocean 
was  continuously  varying  as  the  time  integrations  progressed 
A  100  day  time  mean  zonal  stress  was  calculated  at  a  time 
six  years  prior  to  the  end  of  the  integrations  and  is  shown 
in  Figure  38.   This  stress  was  used  as  a  constant  force 
acting  on  the  ocean  during  the  ocean-only  phase  of  the 
experiment.   Now  it  can  be  interpreted  as  the  mean  stress 
responsible  for  the  stream  function  field  observed  in 
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Figure  33.   The  curl  of  the  stress  has  maxima  at  about  34N 
and  50N  which  are  the  approximate  positions  of  the  centers 
of  the  subtropical  and  subarctic  gyres, respectively.   The 
positions  of  the  maxima  are  about  5  degrees  further  north 
than  actually  observed  in  nature.   This  is  responsible  for 
the  northward  shift  of  the  gyres  and  their  associated 

boundary  currents.   The  stress  varies  from  1.20  to  -1.15 

-2 
dynes  cm   .   These  are  reasonable  values  lending  justifica- 
tion to  the  stress  formulation  discussed  earlier. 

The  state  to  which  the  combined  model  progressed 
was  acceptable  in  its  resemblance  to  the  real  atmosphere- 
ocean  system.   While  it  was  not  the  intention  of  this  work 
to  develop  an  operational  model  for  forecasting  either 
atmospheric  or  oceanic  parameters,  it  is  comforting  to  see 
a  reasonable  distribution  of  the  various  parameters. 
3 .   The  Ocean-Only  Experiment  (Phase  II) 

The  zonal  mean  stress  (Figure  38)  that  had  been 
averaged  over  the  100  day  period  from  day  26,650  to  26,750 
was  used  as  a  constant  forcing  function  for  the  subsequent 
six  year  period  over  which  the  ocean-only  model  equations 
were  integrated.   The  atmospheric  surface  temperature,  T , , 
was  also  averaged  over  this  same  100  day  period.   Figure 
39  shows  a  plot  of  the  zonal  mean  temperature  as  a  function 
of  latitude.   The  temperature  varied  from  2.7  degrees 
centigrade  in  the  north  to  23.3  degrees  in  the  south.   This 
temperature  was  used  as  an  equivalent  atmospheric  tempera- 
ture to  which  the  sea  surface  temperature  was  compared  to 
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determine  the  vertical  heat  flux.   Except  for  the  fixed 
atmospheric  surface  temperature  which  was  independent  of  x, 
the  same  formulation  was  used  for  the  heating  as  was  used 
in  the  combined  model.   Now,  however,  one  degree  of  freedom 
was  removed  and  the  ocean  in  a  very  short  time  became 
steady,  except  for  a  slow  secular  change. 
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Figure  39.   100  day  time  mean  atmospheric  temperature  (C) 
as  a  function  of  latitude. 


The  instantaneous  fields  of  the  ocean  dependent 
variables  from  day  26,700  were  used  as  initial  conditions 
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for  the  ocean -only  phase  of  the  experiment.   During  the  six- 
year  ocean-only  experiment,  data  fields  for  all  oceanic 
variables  were  saved  every  100  days  and  selected  daily  para- 
meters were  saved  at  the  same  five  points  as  previously 
discussed. 

The  2200  day  (six  year)  integration  of  the  ocean- 
only  equations  proceeded  at  a  rapid  pace.   No  problems  were 
encountered  with  adjustment  to  the  new  constant  forcing 
functions  and  a  quasi-steady  state  was  achieved.   At  the 
start  of  the  2200  day  phase  of  the  experiment  the  volume 
average  temperature  of  the  ocean  was  increasing  at  about 
0.5  degrees  per  century  as  previously  discussed.   Through- 
out the  six  year  period  the  warming  rate  remained  about  the 
same,  with  a  0.5  degree  per  century  warming  still  present 
at  the  end.   This  makes  this  author  believe  that  perhaps 
this  is  as  close  to  an  equilibrium  state  as  could  be 
reached  with  the   model  as  formulated,  without  the  expendi- 
ture of  an  excessive  amount  of  computing  time.   An  examina- 
tion of  the  sea  surface  temperature,  Figure  40,  and  the 
area  average  of  the  temperature  at  that  level  shows  that 
the  sea  surface  temperature  is  not  changing  at  all.   The 
rate  of  change  of  average  temperature  at  the  upper  level 
is  exactly  zero.   However,  the  average  temperature  at  the 
lowest  level  is  changing  at  a  rate  of  0.6  degrees  per 
century,  showing  that  the  diffusive  mechanism  by  which  the 
deeper  water  must  be  changed  is  very  slow. 
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Figure  40.   Sea  surface  temperature  from  ocean-only  model 
at  day  28,900. 
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Figure  40  shows  the  sea  surface  isotherms  at  day 
28,900.   An  examination  of  this  temperature  field  at  day 
28,800  was  made  and  the  two  maps  were  identical.   Therefore, 
it  was  considered  unnecessary  to  derive  time  mean  maps  of 
this  quantity.   The  pattern  is  similar  to  the  sea  surface 
temperature  field  at  the  end  of  phase  I  (Figure  34)  in  that 
a  strong  temperature  gradient  is  present  in  the  western 
boundary  layer  at  about  35N.   The  isotherms  throughout-  the 
central  portion  of  the  domain  are  slightly  more  zonal  than 
in  phase  I  where  there  is  an  x  dependence  in  the  forcing, 
but  the  meridional  gradient  in  this  area  remained  about  the 
same.   The  warm  core  tropical  water  (above  23C)  does  not 
extend  into  the  boundary  layer  but  is  visible  as  a  narrow 
band  of  warm  water  just  north  of  the  southern  boundary. 

The  major  difference  in  the  temperature  fields  be- 
tween the  two  phases  is  in  the  northern  half  of  the  domain. 
The  thermal  ridge  present  in  the  northwest  quadrant  in 
phase  I  is  observed  only  as  a  small  perturbation  in  the  6C 
isotherm  in  this  phase.   The  thermal  trough  to  the  west  of 
this  feature  is  similar  to  that  present  in  the  combined 
model  with  a  region  of  water  with  temperature  between  4  and 
6  degrees  extending  over  a  thousand  kilometers.   The  remain- 
der of  the  isotherms  in  the  north  have  a  more  zonal  appear- 
ance.  There  is  a  small  pocket  of  water  below  2C  in  the 
northernmost  region  in  phase  II  that  is  not  present  in 
phase  I.   The  isotherms  in  the  northeast  sector  vary  signi- 
ficantly between  the  two  phases.   In  phase  I  the  thermal 
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current  flows  into  the  northeast  corner  when  the  strongest 
sinking  motion  occurs.   In  phase  II  the  vertical  shear  cur- 
rent flows  into  the  eastern  wall  on  a  more  eastward  heading. 
The  maximum  sinking  motion  is  located  about  800  to  1000  km 
south  of  its  position  in  phase  I.   This  can  be  verified  by 
comparing  Figure  41,  the  vertical  velocity  at  a  depth  of 
300  meters  in  the  ocean-only  case,  with  Figure  37,  the  same 
quantity  in  the  combined  model.   The  northward  curvature  of 
the  isotherms  in  the  central  portions  of  the  eastern  bound- 
ary is  present  in  both  phases,  but  is  less  pronounced  in 
the  ocean-only  case.   An  examination  of  the  meridional  cur- 
rent in  this  area  showed  that  there  was  weak  northward 
motion  in  the  northeastern  boundary  region  in  both  phases 
which  when  coupled  with  the  downwelling  in  that  area,  was 
responsible  for  the  warm  water  there.   The  temperature 
patterns  for  the  levels  at  500  and  1300  meters  resemble 
their  counterparts  from  phase  I  (Figure  35  and  36)  with  the 
same  basic  variations  as  were  seen  in  the  upper  level.   The 
maximum  in  the  temperature  in  the  northeast  corner  has  moved 
south  about  800km  as  the  region  of   sinking  motion  moved. 
The  isotherms  show  no  unexpected  results  and  will  not  be 
depicted  here  . 

The  vertical  motion  patterns  play  an  important  part 
in  the  adjustment  of  temperature  patterns  in  the  boundary 
regions.   Strong  vertical  motions  occur  there  due  to  the 
thermal  circulation.   A  region  of   sinking  water  near  40N 
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Figure  41.   Vertical  motion  in  cm  day"   at  300  meters  at 
day  28,900  of  phase  II.   Positive  indicates  upward  motion 
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in  the  western  boundary  region  causes  the  weak  gradient  be- 
tween the  12C  and  14C  isotherms.   The  isotherms  along  the 
western  boundary  are  shifted  southward.   That  is,  the  14C 
isotherm  in  the  ocean-only  model  is  at  the  latitude  at 
which  the  intermediate  15C  line  (not  drawn)  was  in  the  com- 
bined model . 

Figure  42  shows  the  stream  function  field  for  the 
ocean-only  model  at  day  28,900.   It  was  also  found  unneces- 
sary to  compute  a  time  mean  of  this  quantity.   The  subarctic 
gyre  and  the  subtropical  gyre  both  appear  to  be  slightly 
further  north  in  the  ocean-only  model  than  they  are  in  the 
combined  model  (Figure  33).   Those  gyres  have  about  the 
same  intensity  in  the  western  boundary  region  as  the  corres- 
ponding gyres  in  phase  I.   This  was  expected  since  a  mean 
stress  from  the  combined  model  was  their  driving  mechanism. 
The  northwestward  flow  on  the  east  side  of  the  subarctic 
gyre  is  still  present  and  causes  some  warm  water  to  be 
advected,  but  at  this  time  the  large  amplitude  ridge  that 
was  observed  in  phase  I  was  not  present.   The  ocean-only 
sea  surface  temperature  appears  to  be  strongly  tied  to  the 
imposed  zonal  atmospheric  temperature,  except  in  the  west- 
ern boundary  region  and  to  a  lesser  degree  in  the  eastern 
boundary  region. 

The  maximum  value  of  the  meridional  current  is  re- 
duced somewhat  in  the  ocean-only  model.   In  the  western 
boundary  region,  the  maximum  meridional  velocity  is  47  cm 
sec    in  the  ocean-only  model  compared  to  51  cm  sec    in 
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Figure  42 
phase  II. 


Stream  function  pattern  at  day  28,900  in 
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the  combined  model.   Since  the  vertical  mean  current  appears 
about  the  same  in  both  models,  the  difference  is  in  the 
vertical  shear  current.   The  difference  in   the  predicted 
ocean  variables  between  the  ocean-only  and  combined  model 
phases  is  small.   Radical  changes  did  not  occur.   The  dif- 
ferences that  are  present  are  due  to  the  removal  of  the 
time  varying  atmospheric  forcing  mechanism  in  both  the  ex- 
change of  momentum  and  heat  in  the  ocean-only  model. 
4 .   The  Atmosphere-Only  Model  (Phase  III) 


The  zonal  part,  T     ,  of  the  combined  model  sea- 
sea 

surface  temperature  was  averaged  over  the  same  100  day 
period  as  were  the  zonal  stress  and  atmospheric  surface 
temperature.   This  quantity  was  used  as  the  underlying  tem- 
perature with  which  the  atmospheric  temperature  was  com- 
pared for  the  calculation  of  the  vertical  flux  of  sensible 

heat.   Figure  43  shows  the  profile  of  the  time  and  longi- 

t 

tudinal  mean  of  sea  surface  temperature,  T     ,  which 

sea 

closely  resembles  the  mean  atmospheric  temperature  shown 
in  Figure  39 . 

The  initial  conditions  for  phase  III  were  the 
Fourier  coefficients,  A  through  F,  from  day  26,700  of 
phase  I.   There  was  essentially  no  adjustment  required  and 
the  2200  day  integration  was  conducted,  with  data  retained 
as  before,  in  a  minimum  of  computing  time.   At  the  end  of 
this  phase  the  mean  state  of  the  atmosphere  over  the  last 
100  days  was  examined  along  with  the  state  of  the  atmosphere 
from  day  28,900.   These  maps  were  compared  with  the  atmos- 
phere at  the  end  of  phase  I. 
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Figure  43.   100  day  mean  zonal  component  of  sea  surface 
temperature  (C) . 


Figure  44  shows  the  100  day  mean  distribution  of  the 
10  0  0  mb  height  deviation  as  a  function  of  latitude.   Since 
the  thermal  forcing  is  now  a  function  of  y  only,  the  mean 
fields  are  purely  zonal  (with  the  exception  of  small  errors 
due  to  averaging  and  truncation).   The  deviation  in  height 
extends  from  18  meters  at  35N  to  -28  meters  at  54N  and  the 
deviation  pattern  is  such  that  a  geostrophic  wind  regime 
of  easterlies,  westerlies,  easterlies  would  exist  (as  it 
does).   Figure  45  shows  the  100  day  mean  distribution  of 
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Figure  44.   Distribution  of  100  day  mean  1000  rab  height 
deviation  as  a  function  of  latitude. 


the  750  rab  height  deviation  as  a  function  of  latitude.  The 
deviation  goes  from  118  meters  at  the  southern  boundary  to 
-120  meters  at  the  northern  boundary.   The  tight  negative 
height  gradient  through  the  central  latitudes  agrees  with 
the  existence  of  a  regime  of  westerlies  in  that  area. 
Figure  46  depicts  the  distribution  of  the  100  day  mean 
distribution  of  the  250  mb  height  deviation  as  a  function 
of  latitude.   The  gradient  that  exists  from  410  meters  in 
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the  south  to-407  meters  in  the  north  also  agrees  with  the 
westerly  regime  at  that  level. 
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Figure  45.   Distribution  of  100  day  mean  750  rab  height 
deviation  as  a  function  of  latitude. 


Instantaneous  maps  at  day  28,900  of  the  distribution 
of  the  height  deviation  show  a  highly  disturbed  atmosphere. 
Figure  47,  the  1000  mb  height  deviation  field  at  day  28,900 
analyzed  at  25-meter  intervals,  shows  a  pattern  that  has 
a  disturbance  amplitude  which  is  about  the  same  as  its 
counterpart  in  the  combined  model  (Figure  28).   Figure  48 
shows  the  distribution  of  the  250  mb  height  deviation 
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analyzed  at  60  meter  intervals.   This  flow  is  quite  similar 
to  that  at  the  same  level  at  the  same  time  from  the  combined 
model  (Figure  29).   Figure  49  depicts  the  distribution  of 
the  750  mb  height  deviation  analyzed  at  30  meter  intervals. 
The  wave-like  structure  here  has  a  greater  amplitude  and 
the  deviation  is  larger  than  the  same  fields  in  the  combined 
model  (Figure  30) .   The  occurrence  of  the  difference  in  the 
amplitude  of  the  patterns  is  not  surprising  since  the  model 
atmosphere  is  continuously  changing.   It  would  be  a  coin- 
cidence if  the  structure  of  both  the  atmosphere-only  and 
combined  model  atmosphere  were  identical  at  any  given  time. 
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Figure  46.   Distribution  of  100  day  mean  250  mb  height 
deviation  as  a  function  of  latitude. 
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Figure  47.   Distribution  of  1000  mb  height  deviation  at 
day  28,900  in  phase  III. 
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Figure  48.   Distribution  of  250  mb  height  deviation  at 
day  28,900  in  phase  III. 


150 


Longitude 


Figure  49.   Distribution  of  750  rab  height  deviation  at 
day  28,900  in  phase  III. 
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The  100  day  average  zonal    wind  profile  in  the  atmos- 
phere-only model  is  very  similar  to  that  in  the  combined 
model  shown  in  Figure  27.   A  very  slight  reduction  in  the 
jet  maximum  at  250  mb  is  the  only  difference  so  the  atmos- 
phere-only wind  profile  will  not  be  presented  graphically. 

The  time  averaged  results  of  the  atmosphere-only 
phase  of  the  experiment  were  very  similar  to  the  time 
averaged  atmospheric  state  arrived  at  by  the  combined  model 
Since  the  sea  surface  temperature  is  a  slowly  varying  para- 
meter compared  to  the  atmospheric  time  scale,  holding  it 
constant  approximates  long  term  actual  conditions  fairly 
well.   It  would  be  surprising  if  a  radical  change  occurred 
in  the  time  averaged  fields  in  the  atmosphere-only  formu- 
lation. 

5 .   Energy  balance 

An  examination  of  the  energy  in  the  atmosphere  was 
made  at  the  completion  of  phase  I  and  phase  III.   Rather 
than  study  instantaneous  values  of  the  energies  and  energy 
transformation  as  was  done  in  the  earlier  discussion  of  the 
energy  balance,  values  averaged  over  the  last  100  days  of 
the  experiment  were  examined.   The  energy  flow  convention 
as  established  in  Figure  4  was  used. 

The  mean  energy  exchanges  for  the  combined  model 
and  the  atmosphere-only  model  are  shown  in  Figures  50  and 
51,  respectively.   It  is  interesting  to  note  the  similarity 
in  the  flow  diagrams  in  the  two  phases.   The  diabatic  heat- 
ing generates  about  the  same  amount  of  zonal  potential 


152 


Figure  50.   Flow  diagram  of  mean  atmospheric  energy  from 
the  combined  model.   Energy  and  energy  transformation 


units  are  in  watt  m    (  =  10   erg  cm    sec   )  . 


-2 
energy  in  both  cases,  1.16  watt  m    in  phase  I  and  1.19 

-2 
watt  m    in  phase  III.   The  conversion  of  zonal  potential 

to  eddy  potential  and  then  into  eddy  kinetic  energy  is 

-2 
identical  in  both  phases  at  a  value  of  1.26  watt  m   .   The 

conversion  of  eddy  kinetic  to  zonal  kinetic  energy  is  also 
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Figure  51.   Flow  diagram  of  mean  atmospheric  energy  from 
the  atmospheric-only  model.   Same  units  as  in  Figure  50. 


very  similar  in  the  two  phases  with  values  of  0.37  and  0.33 
watt  m   .   The  dissipation  of  the  kinetic  energy  (both 
zonal  and  eddy)  by  both  internal  and  surface  friction  are 
also  very  similar.   The  only  conversion  that  differs  by  a 
significant  percentage  is  the  conversion  of  zonal  kinetic 
to  zonal  potential  energy.   However,  this  has  been  pointed 
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out  earlier  as  being  a  very  small  quantity  and  no  comment 
on  the  difference  between  the  two  phases  is  offered  here. 

It  is  not  surprising  that  these  flow  diagrams  re- 
semble each  other.   The  only  difference  in  the  two  models 
lies  in  the  method  by  which  heat  is  exchanged  between  the 
atmosphere  and  ocean  or  underlying  wet  surface  as  the  case 
may  be.   In  the  calculation  of  the  sensible  heat  flux 
necessary  for  the  generation  of  zonal  potential  energy  in 
the  combined  model  only  the  zonal  part  of  the  sea  surface 
temperature  was  used.   Therefore,  when  the  disturbance 
part  of  the  sea  surface  temperature  was  neglected  in  the 
atmosphere-only  model,  it  did  not  affect  the  generation  of 
energy  very  much.   In  the  other  energy  exchanges  the  differ- 
ences in  the  disturbance  quantities  were  not  sufficient  to 
make  much  difference. 

These  diagrams  compare  favorably  with  Oort's  dia- 
gram (Figure  6)  as  did  the  earlier  energy  diagram  (Figure 
5).   The  dissipation  of  eddy  kinetic  energy  is  still  too  low 
and  no  new  explanation  is  offered  for  this.   It  is  somewhat 
surprising  that  the  instantaneous  values  of  the  energy  ex- 
change shown  in  Figure  5  compare  so  favorably  with  Figures 
50  and  51  in  so  many  ways.   The  fact  that  the  time  rate  of 
change  of  kinetic  and  potential  energy  (both  zonal  and 
eddy)  were  other  than  zero  in  that  case  shows  that  either 
there  were  small  errors  in  calculations,  or  an  exact  steady 
state  in  the  energies  did  not  exist,  as  was  stated  in  the 
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earlier  discussion.   This  is  in  accord  with  present  know- 
ledge of  the  real  atmosphere  which  is  never  in  exact  steady 
state  . 

B.   SPECTRUM  ANALYSIS 

This  section  will  briefly  discuss  the  spectrum  analysis 
techniques  used  in  this  study.   A  computer  program  in  the 
BIOMEDICAL  series  (BIMED  02T)  developed  at  the  UCLA  Medical 
Center  (Dixon  1970)  was  modified  and  utilized  for  all 
spectral  analyses  conducted  on  the  data  generated  in  this 
study  . 

1 .   Basic  Definitions 

Prior  to  discussing  the  procedures  used  or  the 
results  obtained,  a  brief  review  of  the  definition  of  some 
of  the  basic  statistical  parameters  is  given. 

a.  Variance  -  the  mean  square  deviation  from  the 
mean  of  a  variate. 

b.  Standard  deviation  -  the  positive  square  root 
of  the  variance,  used  to  measure  the  degree  of  variability 
of  a  variate . 

c.  Auto-correlation  coefficients  -  linear  correla- 
tion coefficients  between  a  time  series  and  the  same  time 
series  an  interval  of  time  later. 

d.  Auto-covariance  ~    an  auto-correlation  function 
that  has  not  been  normalized  by  the  variance. 

e.  Cross-correlation  -  correlation  between  two 
time  series,  one  of  which  is  at  an  interval  of  time  later. 


156 


f  .   Cross-covariance  -  a  cross-correlation  function 
that  has  not  been  normalized  by  the  variance. 

g.  Power  spectrum  of  time  series  -  a  function  of 
frequency  which  gives  the  contribution  of  oscillations  at 
each  frequency  to  the  variance  of  a  time  series. 

h.   Cospectrum  -  a  function  of  frequency  which 
gives  the  contribution  of  oscillations  at  that  frequency 
to  the  total  cross-covariance  of  two  time  series. 

i.   Quadrature  spectrum  -  a  function  of  frequency 
which  gives  the  contribution  of  oscillations  at  that  fre- 
quency to  the  total  covariance  of  two  time  series  X(t)  and 
Y(t)  obtained  when  all  the  harmonics  of  time  series  X  were 
delayed  by  a  quarter  period  but  the  Y  series  remained  un- 
changed . 

j.   Phase  -  the  relative  phase  of  the  harmonics 
X(t)  and  Y(t)  is  measured  by  the  ratio  of  the  quadrature 
spectrum  to  the  cospectrum. 
2 .   Spectral  Techniques 

Many  investigators  (for  example,  Yanai  et  al,  1968; 
and  Wallace  1971)  have  utilized  spectral  techniques  to 
study  tropospheric  phenomena  primarily  in  areas  of  sparse 
data  where  a  single  station's  time  varying  record  is 
extremely  important.   The  same  type  of  spectral  techniques 
can  be  used  to  study  the  atmospheric  and  oceanic  data  ob- 
tained from  the   numerical  integrations  conducted  in  this 
study.   At  five  selected  points  (Figure  20)  six  years  of 
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data  were  recorded  at  1.125  day  intervals.   At  each  point 
1960  observations  were  obtained.   These  data  were  analyzed 
and  a  power  spectrum  was  obtained  for  most  variables  at  all 
five  points.   Cr os s -spectrum  analyses  were  conducted  on  sev- 
eral variables  yielding  cospectra,  quadrature  spectra  and 
phase  relationships. 

More  attention  was  focused  on  point  C  because  of 
its  position  in  the  western  boundary  region.  It  was  located 
in  an  area  where  meridional  velocities  were  strong,  tempera- 
ture gradients  were  tight  and  upwelling  was  important. 
However,  the  other  points  in  Figure  20  were  not  neglected 
and  comparisons  of  the  resulting  spectral  statistics  between 
points  were  made. 

The  data  were  recorded  at  an  interval,  (A  t)  ,,  of 

d 

1.125  days  (three  oceanic  time  steps)  during  the  three  six- 
year  experiments.   That  is,  the  data  were  available  for 
analysis  from  the  combined,  ocean-only  and  atmosphere-only 
phases  of  the  experiment.   The  decision  to  record  data  at 
1.125  days  was  a  critical  one.   It  would  require  an  exces- 
sive amount  of  storage  space  to  keep  a  longer  record  of 
data  at  shorter  recording  intervals  and  problems  could  be 
encountered  if  the  interval  was  too  long.   It  is  important 
to  ensure  that  the  series  which  is  being  examined  does  not 
contain  substantial  amounts  of  variance  from  frequencies 
outside  the  range  that  can  be  resolved  with  the  spectra. 
Spectral  information  for  frequencies  less  than  one  cycle/ 


2  (  Z)  t)  ,  the  Nyquist  folding  frequency,  can  be  obtained  fro 


in 
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the  data.   Frequencies  higher  than  this  are  "folded"  or 
"aliased"  into  the  lower  frequency  region.   Wallace  (1971) 
points  out  that  the  high  frequency  limit  does  not  appear  to 
pose  any  serious  problem  from  aliasing  because  there  is 
little  variance  associated  with  frequencies  higher  than  0.5 
cycles  per  day  in  the  intermediate  to  large  scale  disturb- 
ances.  From  experimental  work  Wallace  concluded  that  there 
were  only  minor  differences  in  spectra  computed  with  daily 
or  twice-daily  data.   Since  the  data  examined  in  this  study 
is  of  the  daily  variety,  his  conclusions,  along  with  the 
fact  that  higher  frequencies  due  to  advection  of  short 
waves  cannot  exist  in  the  model  atmosphere  leads  the  author 
to  believe  that  the  time  interval  chosen  for  this  data  re- 
cord was  a  sound  one. 

In  the  BIMED  02T  program,  the  power  spectrum  is  ob- 
tained by  performing  a  harmonic  analysis  on  the  auto-co- 
variance  function  of  the  quantity  to  be  studied.   After 
each  spectral  estimate  was  multiplied  by  the  interval  of 
frequency  (A  f)  of  the  band  width,  the  spectral  estimates 
for  that  frequency  band  can  be  summed  to  obtain  the  contri- 
bution of  the  fluctuations  within  that  specific  frequency 
band  width  to  the  variance  of  the  time  series  which  is 
under  examination. 

3 .   Results  of  Spectral  Investigation 

The  primary  point  of  interest  in  the  spectral  in- 
vestigation was  point  C  (Figure  20)  in  the  western  bound- 
ary region.   A  thorough  investigation  of  atmospheric  and 


159 


oceanic  variables  was  conducted  at  this  point  with  the  fol- 
lowing results . 

a.   Study  of  data  at  point  C. 

The  data  obtained  from  phase  I  (combined  model) 
of  the  experiment  were  examined  first.   The  Fourier  coeffi- 
cients from  the  atmospheric  part  of  the  combined  model  were 
analyzed  and  Figure  52  shows  the  power  spectrum  of  A  and  B 
(solid  line)  and  C  and  D  (dashed  line),  obtained  from  Equa- 
tions (2.22)  to  (2.25),  plotted  on  a  linear  scale  against 
frequency  (cycles  per  day) .   The  power  spectra  of  A  and  B 
are  identical  as  are  the  spectra  of  C  and  D.   The  amplitude 
of  the  power  spectra  of  C  and  D  is  slightly  less  than  10% 
of  that  of  A  and  B.   There  is  a  pronounced  peak  at  about 
7.5  days  in  the  power  spectra  for  these  quantities.   This 
was  expected  since  the  phase  speed  of  the  imposed  baroclinic 
disturbance  in  the  combined  model  was  estimated  at  600  km 
per  day,  a  phase  speed  that  would  yield  a  period  of  8  days. 
In  all  likelihood,  the  estimate  of  the  phase  speed  of  the 
disturbance  was  slightly  low.   The  7.5  day  period  indicates 
that  640  km/day  was  the  speed  with  which  the  disturbance 
moved  across  the  4800  km  wide  domain. 

Figure  53  shows  the  power  spectrum  of  atmos- 
pheric disturbance  kinetic  energy.   There  are  some  small 
peaks  in  the  spectra  in  the  frequency  band  from  .08  to  .24 
cycles  per  day   (12.5  to  4.2  days)  but  the  major  peaks  are 
clearly  at  the  lower  frequencies.   There  are  peaks  at 
periods  of  about  32  and  62  days  indicating  cycles  in  the 
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Figure  52.   Power  spectra  of  A,B,C  and  D  at  point  C  in 
phase  I.   Period  of  spectral  peak  indicated  in  days. 
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Figure  53.   Power  spectrum  of  atmospheric  disturbance 
kinetic  energy  in  phase  I.   Period  of  spectral  peaks 
indicated  in  days. 
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atmospheric  energy  on  the  order  of  one  and  two  months.   If 
seasonal  variations  had  been  included  in  the  specification 
of  the  solar  radiation  and  latent  heat  release,  different 
frequencies  of  oscillation  would  probably  be  observed.   The 
annual  fluctuation  in  the  above  mentioned  parameters  would 
probably  be  reflected  in  the  energy  spectrum. 

Figure  54  shows  the  power  pectrum  of  E,  the  zo- 
nal mean  part  of  \p       that  is  proportional  to  the  mean  height 
and  related  to  the  mean  geostrophic  wind  in  the  atmosphere. 
The  quantity  has  a  large  spectral  peak  at  a  period  of  32 
days  and  smaller  peaks  at  25  and  62  days.  An  examination 
of  the  cospectrum  and  phase  between  E  and  the  disturbance 
energy  shows  that  the  disturbance  energy  rise  at  the  32 
day  period  lags  the  increase  in  the  spectrum  of  E  at  that 
period.   This  suggests  that  the  mean  flow  increased  first 
and  because  of  barotropic  amplification  transferred  energy 
to  the  disturbance. 

Figure  55  depicts  the  power  spectrum  of  the 
oceanic  stream  function.   Peaks  in  the  spectra  are  seen  at 
periods  of  7.5  days  (small),  and  32  and  62  days  (large). 
This  is  the  first  quantitative  indication  of  an  air-sea 
interaction  phenomenon  in  the  model.   The  7.5  day  peak  in 
the  power  spectrum  of  the  stream  function  is  caused  by  the 
contribution  of  the  A  and  B  coefficients  to  the  disturbance 
component  of  the  atmospheric  wind  stress  (Figure  52). 
Since  the  disturbance  components  of  the  stress  are  small  the 
fluctuations  they  cause  are  small  and  therefore  the  peak  in 
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Figure  54.   Power  spectrum  of  E  at  point  C  in  phase  I 
Period  of  spectral  peaks  indicated  in  days. 
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Figure  55.   Power  spectrum  of  ocean  stream  function  at 
point  C  in  phase  I.   Period  of  spectral  peaks  indicated  in 
days  . 
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the  spectrum  is  smaller  at  this  frequency.   The  larger  peaks 
at  about  32  and  62  days  are  caused  by  the  oscillations  in 
the  zonal  component  of  the  wind  stress  which  is  proportional 
to  E.   An  examination  of  the  cospectrum  and  phase  between  E 
and  the  ocean  stream  function  shows  that  the  increase  in  the 
mean  atmospheric  flow  occurs  3.5  days  prior  to  the  fluctua- 
tion in  the  stream  function  at  the  32  day  period.   Since  E 
had  spectral  peaks  at  these  frequencies  it  is  not  surprising 
that  peaks  in  the  same  region  of  the  spectrum  appear  in  the 
ocean  stream  function.   Figure  56  shows  a  100  day  time 
series  of  E  and  the  oceanic  stream  function  obtained  by 
averaging  the  recorded  data  over  seven  day  intervals.   The 
averaging  was  done  to  remove  the  7.5  day  oscillation  from 
the  data  to  examine  the  longer  period  fluctuations  in  more 
detail.   It  is  easy  to  see  that  the  long  period  oscillations 
are  present  in  both  E  and  the  stream  function  and  that  the 
maxima  in  E  leads  the  maxima  in  the  stream  function  by 
several  day s . 

Figure  57  shows  the  power  spectrum  of  the  zonal 
mean  of  the  atmospheric  surface  temperature,  T, .   Since  T, 
is  proportional  to  F,  the  peaks  in  the  spectrum  are  also 
representative  of  peaks  in  the  spectrum  of  F.   This  was 
verified  by  spectrally  analyzing  F  separately,  the  plot  of 
which  is  not  presented.   There  are  peaks  in  T,  of  almost 
equal  amplitude  at  7.5  and  32  days,  and  smaller  peaks  at 
periods  of  25,  62  and  111  days.   Figure  58  shows  the  power 

spectrum  of  T    ,  the  sea  surface  temperature.   Small  peaks 
sea  r 
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Figure  57.   Power  spectrum  of  atmospheric  surface  tempera' 
ture  at  point  C  in  phase  I.   Period  of  spectral  peaks 
indicated  in  days. 
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Figure  58.   Power  spectrum  of  sea  surface  temperature  at 
point  C  in  phase  I.   Period  of  spectral  peaks  indicated 
in  day s . 
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can  be  seen  at  periods  of  about  25  and  32  days  and  a  larger 
peak  at  a  period  of  about  62  days.   These  peaks  coincide 
with  some  of  the  peaks  in  the  spectrum  of  T , .   Figure  59 
shows  the  power  spectrum  of  the  vertical  heat  flux  between 
the  atmosphere  and  the  ocean.   Peaks  in  the  spectrum  occur 
at  7.5,  32  and  62  days.   The  longer  period  peaks  are  at  the 
same  period  as  the  peaks  in  the  atmosphere  and  ocean  surface 
temperatures,  indicating  an  exchange  of  heat  due  to  a  tem- 
perature difference  between  the  two  media.   A  comparison  of 
the  time  series  of  the  atmospheric  surface  temperature,  the 
vertical  heat  flux  and  the  sea  surface  temperature  obtained 
by  averaging  over  seven  day  periods  is  shown  in  Figure  60. 
It  is  readily  observed  that  oscillations  on  the  order  of 
30  days  are  present  in  all  three  series  and  not  so  obvious 
fluctuations  on  the  order  of  60  days  are  present  in  the  sea 
surface  temperature  and  the  vertical  heat  flux.   It  can  be 
seen  that  peaks  in  the  sea  surface  temperature  led  peaks 
in  the  vertical  heat  flux  which  in  turn  led  peaks  in  the 
atmospheric  surface  temperature. 

An  examination  of  the  cospectrum  and  phase  be- 
tween E  and  T     shows  that  the  increase  in  E  at  about  32 
sea 

days  precedes  the  fluctuation  in  sea  surface  temperature  by 

11  days.   This  was  perhaps  a  result  of  an  increase  in  the 

mean  conal  wind  stress  causing  an  increase  in  the  wind 

driven  part  of  the  ocean  circulation  and  an  increase  in  the 

advection  of  heat  from  the  south  at  point  C.   An  examination 

of  the  cospectrum  of  T     and  the  ocean  stream  function 

sea 
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Figure  59.   Power  spectrum  of  vertical  heat  flux  at  point  C 
in  phase  1.   Period  of  spectral  peak  indicated  in  days. 
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shows  that  the  fluctuations   in  the  stream  function  at  32 
days  occurs  7.6  days  prior  to  the  oscillation  in  T 

5  6a 

These  phase  relationships  coupled  with  the  previously  stated 
facts  leads  to  the  following  conclusions.   The  oscillations 
which  were  reflected  as  peaks  in  the  power  spectra  of  the 
variables  were  caused  by  an  air-sea  interaction  phenomenon 
and  were  not  merely  spurious  reflections  of  energy  due  to 
sampling  or  calculating  problems.   An  explanation  of  the 
probable  series  of  events  causing  these  interactions  is  as 
follows.   The  increase  in  the  mean  atmospheric  flow,  re- 
flected in  E,  caused  an  increase  in  the  wind  driven  part 
of  the  ocean  circulation  3.5  days  later.   This  increased 
circulation  brought  about  an  increase  in  the  advection  of 
heat  from  the  south  at  point  C  and  an  oscillation  in  the 
sea  surface  temperature,  at  a  period  of  32  days,  occurred 
7.6  days  later.   The  total  time  from  the  increase  in  E  to 

the  increase  in  T     was  about  11  days.   The  increase  in 

sea 

T     caused  an  increase  in  the  vertical  heat  flux  and  an 
sea 

increase  in  the  atmospheric  surface  temperature.   An  exam- 
ination of  the  cospectrum  and  phase  of  T     and  T,  shows 

sea       4 

that  the  oscillation   in  T     at  32  days  leads  the  rise  in 

sea  J 

the  spectrum  of  T,  by  6.2  days.   An  exchange  of  energy  be- 
tween T  (F)  and  E  is  then  required  to  maintain  the  cycle. 
4 

An  examination  of  the  cospectrum  and  phase  of  E  and  F  shows 
that  the  peak  in  F  occurs  prior  to  the  peak  in  E  at  32  days. 
When  these  interactions  are  taking  place  in  the  combined 
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model  and  oscillations  at  various  frequencies  are  present, 
almost  every  dependent  variable  in  the  model  will  be  affec- 
ted in  some  way.   It  has  been  shown  that  the  long  period 
oscillations  which  are  present  in  the  coupled  model  include 
the  ocean  as  an  important  part  in  the  oscillation  cycle. 
It  should  not  be  surprising  to  find  that  these  long  period 
oscillations  do  not  exist  when  either  the  atmosphere  or 
ocean  is  held  constant  in  time. 

The  data  from  phase  II  of  the  experiment  were 
examined  next.   With  the  time  varying  forcing  functions 
removed,  the  ocean  became  steady  in  a  short  time  with  no 
significant  peaks  in  any  of  the  spectra.   Figure  61  shows 

the  power  spectrum  for  T    .   This  is  typical  of  all  the 

sea  J 

ocean  variables  in  this  phase  in  that  the  power  in  the 
spectrum  is  in  the  lowest  frequencies.   A  comparison  of 
Figure  61  and  58  provides  a  convincing  argument  that  the 
time  varying  atmospheric  circulation  is  responsible  for  the 
transient  responses  of  the  ocean  surface  temperature  at 
periods  of  about  one  and  two  months. 

An  interesting,  and  initially  perplexing,   peak 
appeared  in  the  ocean  stream  function  at  7.5  days  in  the 
ocean-only  data  (not  shown).   The  ocean-only  model  required 
a  period  of  time  to  adjust  to  the  new  forcing  and  during 
that  time  the  stream  function  oscillated  and  then  became 
constant.  Although  the  energy  in  the  spectrum  was  diluted 
by  the  long  period  of  constancy,  the  initial  oscillation  made 
itself  evident  at  the  frequency  at  which  it  fluctuated. 
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Figure    61.       Power    spectrum    of  sea    surface    temperature    at 
point    C    in    phase    II. 
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The  data  from  phase  III  (atmosphere-only)  were 
analyzed  and  Figure  62  shows  the  power  spectra  of  A  and  B 
(solid  line)  and  C  and  D  (dashed  line) .   The  four  quantities 
all  have  a  peak  at  a  period  of  8.4  days.   This  is  different 
than  the  period  at  which  there  is  a  peak  in  the  spectra  of 
the  same  quantities  in  the  combined  model.   For  the  4800  km 
wavelength  used  in  this  study,  a  period  of  8.4  days  yields 
a  phase  velocity  of  570  km  day   .   This  is  about  12%  slower 
than  the  phase  speed  for  the  atmospheric  wave  in  the  com- 
bined model.   By  eliminating  the  zonal  variation  in  the  sea 
surface  temperature,  the  sensible  heat  exchange  was  changed. 
A  comparison  of  the  combined  model  and  atmosphere-only  model 
shows  that  the  thermal  wind  was  reduced  in  the  latter. 
This  reduction  in  the  thermal  wind,  shown  in  Table  4,  was 
influential  in  retarding  the  phase  speed  of  the  disturbance. 

The  power  spectrum  of  the  disturbance  kinetic 
energy  is  shown  in  Figure  63.   A  large  peak  at  the  6.7  per- 
iod is  observed  along  with  a  smaller,  wider  peak  in  the  45 
to  50  day  period  range.   The  peak  at  the  6.7  day  period, 
which  is  near  the  8.4  day  period  in  the  disturbance  quanti- 
ties, is  perplexing  in  light  of  the  fact  that  there  is  no 
sharp  peak  around  the  7.5  day  peak  in  the  disturbance  kine- 
tic energy  (Figure  53)  in  the  combined  model.   A  possible 
explanation  is  that  the  6.7  day  peak  is  very  sharp  and 
while  it  has  a  considerable  amplitude,  the  total  contribu- 
tion to  the  variance  of  the  disturbance  is  not  much  greater 
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Figure  62.   Power  spectrum  of  A,B,C  and  D  at  point  C  in 
phase  III.   Period  of  spectral  peaks  indicated  in  days. 
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LATITUDE 

THERMAL  WIND 

THERMAL  WIND 

(I) 

(III) 

72.0 

0.0 

0.0 

68.4 

1.2 

3.0 

64  .8 

4.2 

3.7 

61.2 

3.9 

3.5 

57.6 

4.6 

4.3 

54  .0 

5.7 

5.6 

50.4 

6.6 

7  .5 

46.8 

8.5 

9.4 

43.2 

10.2 

9.2 

39.2 

9.1 

7.5 

36  .0 

6.8 

5.6 

32.4 

4.4 

3.9 

28.8 

3.0 

3.1 

25  .2 

2.3 

2.8 

21.6 

1.3 

2.5 

18.0 

0.0 

0.0 

TABLE  4.   Comparison  values  of  mean  atmospheric  thermal 
wind  from  combined  model  (I)  and  atmosphere-only  model 
(III),  (in  m  sec"l). 
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Figure  63.   Power  spectrum  of  atmospheric  disturbance 
energy  in  phase  III. 
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than  the  small  diffuse  indication  of  a  peak  around  the  7.5 
day  period  in  the  combined  case.   In  the  spectra  for  the 
quantities  E  and  F,  peaks  at  these  same  periods  are  ob- 
served, showing  the  exchange  occurring  between  the  disturb- 
ance and  mean  flow. 

An  examination  of  data  in  phases  I,  II  and  III 
at  point  C  indicates  that  the  peaks  in  the  spectra  noted  at 
about  one  and  two  months  in  the  combined  model  depend  on 
atmosphere-ocean  interactions  for  their  existence  and  are 
not  simply  spurious  or  caused  by  aliasing.   The  fact  that 
these  interactions  are  on  the  order  of  one  and  two  months 
is  not  the  most  significant  point,  because  the  period 
would  undoubtedly  change  with  some  other  period  of  atmos- 
pheric forcing.   What  is  significant  is  that  these  longer 
period  fluctuations  can  be  attributed  to  air-sea  interac- 
tions.  The  absence  of  such  peaks  in  the  spectra  in  either 
the  atmosphere-only  or  ocean-only  phases  shows  that  the 
fluctuations  causing  the  peaks  are  phenomena  occurring  only 
when  the  freedom  found  in  the  coupled  model  is  allowed. 
Although  there  is  a  peak  in  both  atmospheric  and  oceanic 
variables  in  the  coupled  model  at  7.5  days,  this  is  not 
considered  as  significant  as  the  peaks  at  the  longer  per- 
iods because  this  frequency  of  oscillation  was  forced  into 
the  model  with  the  selection  of  the  4800  km  wavelength  for 
the  single  wave  in  x.   Changing  this  wavelength  would  cause 
the  peak  in  the  spectrum  at  7.5  days  to  shift  and  undoubted' 
ly  would  have  many  other  effects  on  the  model.   It  would 
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be  interesting  to  see  if  the  peaks  in  the  spectra  at  the 
lower  frequencies  would  be  shifted. 

Spectral  analyses  of  the  kinetic  energy  in  the 
ocean  of  the  combined  model  at  the  three  ocean  levels  were 
made  at  point  C.   An  approximation  to  the  total  kinetic 
energy  at  point  C  was  used  for  this  study.   The  total  kine^ 
tic  energy  should  include  the  vertical  mean  kinetic  energy 
and  the  vertical  shear  kinetic  energy.   A  mathematical  for' 
mulation  for  this  is 


K'E       =  T[ (u' 
ocean    I 


!V  +  (v,  +  91)2] 
3y;     K  dxJ     J 


(6.2) 


In  recording  the  data  at  point  C  no  provisions  were  made 

for  the  calculation  of  -r —   and  so  the  zonal  component  of 

3y 

the  vertical  mean  current  at  this  point  was  neglected.  The 
zonal  component  of  the  vertical  mean  current  was  small  and 
therefore  this  was  not  a  bad  assumption.   The  kinetic 
energy  computed  at  each  oceanic  level  was 


K*E       =  y[(u')2  +  (v»  +  IV]. 
ocean    2  dx 


(6.3) 


Figure  64  shows  the  power  spectrum  of  kinetic 
energy  at  100  meters  in  the  ocean  of  the  combined  model 
versus  frequency  in  cycles  per  day,  both  plotted  on  log 
scales.   Figure  65  is  the  same  plot  for  500  meters  and 
Figure  66  shows  the  1300  meter  depth  analysis.   Peaks  in 
the  spectra  at  periods  of  62,  25-32,  and  7.5  days  are  ob- 
served at  all  three  levels.   This  is  consistent  with  earl- 
ier findings.   The  amplitude  of  the  25-32  day  period 
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Figure  64.  Power  spectrum  of  oceanic  kinetic 
meters  at  point  C  in  phase  I.  Periods  of  spe 
indicated  in  days. 
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Figure  65.   Power  spectrum  of  ocean  kinetic  energy  at  500 
meters  at  point  C  in  phase  I.   Periods  of  spectral  peaks 
indicated  in  days. 
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Figure  66.  Power  spectrum  of  ocean  kinetic  energy  at  1300 
meters  at  point  C  in  phase  I.   Periods  of  spectral  peaks 
indicated  in  days. 
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is  the  largest  at  the  upper  level  while  the  amplitude  of 
the  62  day  period  and  25-32  day  period  are  about  the  same 
at  the  two  lower  levels.   A  comparison  of  Figures  65  and  66 
shows  that  the  spectra  for  the  two  lowest  levels  in  the 
model  nearly  coincide.   The  maximum  amplitude  in  the  power 
spectra  decreased  by  an  order  of  magnitude  between  the  100 
meter  level  and  the  500  meter  level.   This  lessening  in  the 
variance  of  the  kinetic  energy  between  these  levels  is  not 
surprising  since  the  lower  levels  are  less  active.   The 
amplitude  remained  about  the  same  between  the  500  and  1300 
met  er  level s . 

Thompson  (1971)  analyzed  a  long  series  of  data 
recorded  by  the  Woods  Hole  Oceanographic  Institution  at  a 
site  located  at  39  20N,  70W  in  the  North  Atlantic.   In  the 
model  this  would  be  analogous  to  a  position  approximately 
200  km  northeast  of  point  C.   The  data  were  recorded  at 
900  second  intervals.   The  folding  frequency  associated 
with  this  small  recording  interval  allowed  Thompson  to  in- 
vestigate the  shorter  period  oscillations.   Some  gaps  in 
the  data  forced  the  estimation  of  some  of  the  spectra  in 
the  low  frequency  region.   He  analyzed  the  northward  and 
eastward  components  of  velocity  and  through  a  summing  and 
averaging  process  calculated  the  "kinetic  energy  density" 
at  several  depths  in  the  ocean  (not  shown).   None  of  these 
spectra  showed  any  peaks  at  periods  greater  than  one  day, 
although  osciallations  with  periods  of  about  ten  hours  were 
observed.   The  spectra  for  the  three  lowest  depths  were 
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nearly  coincident,  and  there  was  a  large  decrease  in  the 
amplitude  of  the  spectra  between  100  meters  and  500  meters. 
In  this  feature  the  study  of  Thompson  and  the  study  con- 
ducted in  this  work  were  similar. 

b.   Study  of  data  at  point  A. 

At  point  A  (northeast  area)  the  power  spectrum 

of  the  coefficients  A,B,C  and  D  were  greatly  reduced  [a 

2    -12 
maximum  value  of  16  (cm  sec   )   cycles  per  day  for  A  and  B 

as  compared  to  675.0  at  point  C],  although  the  peak  in  the 
spectra  was  at  the  same  period,  7.5  days.   An  examination 
of  the  coefficient,  E,  shows  the  spectral  peak  at  32  and  62 
days  and  a  large  peak  at  120  days.   Looking  back  at  Figure 
53,  a  smaller  peak  in  the  region  of  the  120  day  period  can 
be  seen  also.   The  oceanic  stream  function  at  point  A  shows 
oscillations  that  have  smaller  amplitudes  in  the  power 
spectrum  than  do  those  of  the  same  quantity  at  point  C. 
The  biggest  difference  between  the  power  spectra  of  vari- 
ables at  point  A  and  those  at  point  C  occurred  with  T 

sea 

There  are   no  peaks  at  the  32  or  62  day  periods  and  all  the 
power  is  in  the  lowest  frequency,  much  the  same  as  seen  in 
Figure  58,  the  same  quantity  in  the  ocean-only  case. 

At  point  A  in  the  ocean-only  phase,  the  power 
spectra  resembled  those  from  the  same  phase  at  point  C. 
The  maximum  of  the  power  spectrum  that  existed  in  the  low- 
est frequency  was  increased  but  will  not  be  discussed  fur- 
ther since  it  adds  no  new  information.   An  increase  in  the 
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power  in  the  lowest  frequency  merely  indicates  that  the  long 
terra  change  was  larger.   In  the  atmosphere-only  phase,  the 
power  spectra  for  the  Fourier  coefficients  resemble  those 
at  point  C,  except  that  the  maximum  of  the  spectra  at  the 
8. A  day  period  is  reduced  to  approximately  10%  of  the  value 
it  had  at  point  C.   At  point  A,  the  meridional  temperature 
gradient  is  rather  weak  and  the  amplitude  of  the  baroclinic 
wave  is  much  less  than  in  the  lower  latitudes  of  point  C. 
c.   Study  of  data  at  Point  B,  D  and  E 

At  point  B  (central  area)  the  power  spectra  of 
the  coefficients  A,  B,  C,  and  D  were  slightly  larger  than 
those  at  point  C.   This  location  at  45N  is  in  an  area  of 
large  thermal  wind,  so  the  amplitude  should  be  greater. 
Except  for  the  increase  in  amplitudes  the  spectra  are  very 
similar  to  those  of  the  variables  at  point  C  (Figures  52 
through  56) .   The  power  spectrum  of  T    ,  though,  is  more 

5  6  ct 

like  that  at  point  A,  with  all  the  power  at  the  lowest  fre- 
quency. This  indicates  extremely  long  period  fluctuations, 
probably  too  long  for  this  data  record  to  portray. 

The  ocean-only  variables  at  point  B  resemble 
those  at  point  A  and  will  not  be  discussed.   The  Fourier 
coefficients,  A  through  D,  have  a  peak  at  the  same  8.4  day 
period  and  the  maximum  in  the  power  spectrum  is  increased. 
No  new  information  was  obtained  from  point  B. 

A  separate  discussion  of  the  data  at  point  D 
will  not  be  made  because  of  similarity  of  results  to  point 
B.   Strange  results  were  obtained  when  the  data  at  point  E 
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were  spectrally  analyzed.  The  Fourier  coefficients,  A 
through  D,  had  multiple  peaks  in  the  spectra  that  were  not 
representative.  It  was  decided  that  this  point  was  selec- 
ted too  close  to  the  southern  boundary  and  the  artificial- 
ity of  the  wall  caused  unusual  oscillations  in  the 
atmospheric  data.  The  results  from  point  E  were  discarded 
with  no  further  comment. 

The  results  of  the  spectral  study  were  that  the 
ocean  circulation  caused  transients  to  appear  in  the  atmos- 
phereic  circulation  in  the  one  and  two  month  periods  and 
possibly  with  a  period  of  120  days  and  that  similar  tran- 
sients appeared  in  the  ocean  circulation  due  to  the  pre- 
sence of  the  fluctuating  atmosphere.   A  great  deal  of 
information  has  been  wrested  from  the  spectral  statistics 
discussed  here,  but  a. great  deal  more  could  be  obtained  if 
time  permi  tted  . 


VII .   SUMMARY 

A  numerical  model  of  a  combined  atmospheric-ocean  system 
was  developed  for  use  in  an  investigation  of  the  long-term, 
large-scale  air-sea  interactions  and  the  effect  of  these  in- 
teractions on  transients  that  appear  in  the  atmosphere  and 
oceanic  circulations.   The  problems  associated  with  observa- 
tional studies  of  the  air-sea  system  were  eliminated  by 
using  a  numerical  model  that  could  produce  a  long-term  simu- 
lation of  an  interacting  system  in  a  reasonable  amount  of 
computing  time. 

The  atmospheric  model  used  in  this  study  was  fashioned 
after  the  two  level,  quasi-geos trophic  model  developed  by 
Phillips  (1956).   New  dependent  variables  were  introduced 
by  the  use  of  a  one-dimens ional  (in  x)  Fourier  series,  in 
which  a  single  mode   was  retained,  to  represent  the  origi- 
nal variable.   This  method  of  approach  allowed  the  rapid 
integration  of  the  atmospheric  equations  since  the  number 
of  degrees  of  freedom  were  greatly  reduced.   Without  this 
modification  the  combined  model  could  not  have  been  inte- 
grated for  a  sufficiently  long  time  period  to  obtain  mean- 
ingful results  . 

The  ocean  model  utilized  in  this  work  was  a  modified 
version  of  Haney's  (1971a)  primitive  equation  model.   The 
modifications  included  a  reduction  in  the  number  of  levels, 
a  reduction  in  horizontal  region  size  (and  the  number  of 
grid  points),  a  reduction  in  the  depth  of  the  basin  and  the 
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conversion  to  3~plane  (cartesian)  geometry.   The  primary 
reason  for  these  modifications  was  to  increase  the  ratio  of 
simulation  time  to  computing  time  and  permit  a  long  term 
simulation  which  was  accomplished. 

The  coupling  of  the  model  was  accomplished  through  the 
exchange  of  momentum  through  the  wind  stress  and  the  ex- 
change of  heat  through  the  turbulent  exchange  process  which 
was  proportional  to  the  air-sea  temperature  difference. 
Assumptions  were  made  to  simplify  both  of  these  exchanges. 
The  atmospheric  model,  however,  had  to  be  restricted  to 
simplified  formulations  if  it  was  to  be  integrated  success- 
fully for  the  necessary  length  of  time.   The  absence  of 
moisture  in  the  atmospheric  model  required  further  simpli- 
fication in  the  formulation  of  the  latent  heat  flux.   These 
simplifications  were  indeed  necessary  and  an  examination  of 
the  results  showed  that  the  model  produced  atmospheric  and 
oceanic  circulations  that  were  adequate  for  this  experiment 
in  their  portrayal  of  the  real  atmosphere-ocean  system. 

The  experiment  was  conducted  in  three  phases.   Phase  I 
was  a  long  term  integration  of  the  combined  model.   The 
computing  time  was  reduced  by  the  use  of  a  larger  grid  size 
for  the  first  portion  of  this  phase.   The  400  km  grid  dis- 
tance was  that  used  to  accelerate  the  integration  of  the 
combined  model  equations  in  the  first  part  of  Phase  I  was 
just  adequate  to  resolve  the  western  boundary  current. 
The  coefficient  of  lateral  diffusion  of  heat  had  been  re- 
duced to  prevent  large  horizontal  transports  of  heat  by- 


190 


diffusion  and  keep  the  meridional  temperature  gradient  at  a 
realisticly  high  level.   The  lateral  diffusion  of  momentum 
was  reduced,  but  could  not  be  made  small  enough  so  that  a 
truly  nonlinear  solution  could  be  obtained  without  having 
the  width  of  the  boundary  current  become  smaller  than  the 
horizontal  grid  distance.   The  unrealistic  spatial  oscilla- 
tions seen  in  the  results  obtained  with  the  400  km  grid  were 
a  direct  result  of  the  large  grid  size. 

The  spatial  oscillations  observed  in  the  400  km  grid 
model  results  quickly  disappeared  in  the  200  km  grid  version 
of  the  combined  model.   A  state  of  quasi-equilibrium  was 
reached  in  the  long  term  integrations  of  the  combined  model. 
The  results  shown  throughout  this  thesis  were  considered  to 
be  satisfactory  for  the  purpose  of  investigating  the  effect 
of  the  oceanic  (atmospheric)  circulation  on  the  time-depend- 
ent behavior  of  the  atmosphere  (ocean) .   Some  features  in 
the  results  considered  worthy  of  note  are  summarized  below. 

A  western  boundary  current  developed  with  very  strong 
northward  flow,  analogous  to  the  Gulf  Stream,  on  the  west 
side  of  the  sub-tropical  gyre  and  strong  southward  flow, 
analogous  to  the  Labrador  Current,  on  the  west  side  of  the 
sub-arctic  gyre.   The  gyres  themselves  are  positioned  as 
found  in  nature  and  are  a  result  of  the  surface  wind  stress 
developed  in  the  atmospheric  model.   The  conditions  in  the 
northeast  corner  of  the  domain  are  similar  to  those  found 
in  the  northeast  Pacific.   The  warmer  water  is  due  to  the 
convergence  of  the  vertical  shear  current  into  the 
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northeastern  boundary  resulting  in  downwelling  in  that  re- 
gion.  The  region  of  warm  water  and  weak  currents  in  the 
south  central  region  of  the  domain  is  analogous  to  the  water 
found  in  the  Sargasso  Sea  area  of  the  Atlantic. 

The  zonal  winds  in  the  atmosphere  develop  a  horizontal 
and  vertical  structure  similar  to  that  found  in  the  real  at- 
mosphere.  A  three  cell  meridional  circulation  was  formed 
and  slowly  meandered  several  degrees  of  latitude  north  and 
south  of  its  mean  position. 

The  second  phase  of  the  experiment  consisted   of  a  six 
year  integration  of  the  ocean-only  model  and  a  state  of 
equilibrium  was  quickly  reached.   The  upper  level  of  the 
ocean  had  reached  a  state  of  almost  steady  temperature  but 
the  lower  level  still  showed  a  very  slow  increase  in  the 
temperature.   This  was  a  good  example  of  the  extremely  long 
time  necessary  for  the  vertical  diffusion  mechanism  to  bring 
the  deep  water  into  thermal  equilibrium. 

Phase  III  of  the  experiment  consisted  of  the  integration 
of  the  atmosphere-only  model  equations.   Energy  balance 
studies  of  the  atmosphere  from  this  phase  were  made  and  com- 
pared with  the  same  studies  from  phase  I.   The  energy 
studies  compared  favorably  with  each  other  and  with  the 
existing  knowledge  of  atmospheric  energy  exchanges. 

Data  produced  by  the  combined  model  and  recorded  at 
selected  points  were  spectrally  analyzed  yielding  peaks  in 
the  spectra  of  both  oceanic  and  atmospheric  variables  at 
periods  on  the  order  of  one  and  two  months.   Cross-spectrum 
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analyses  and  time  series  analyses  of  several  atmospheric  and 
oceanic  variables  led  to  the  following  conclusion  about  the 
air-sea  interactions.   The  longer  period  oscillations  were 
present  in  the  data  and  were  attributed  to  air-sea  interac- 
tion.  An  increase  in  the  mean  atmospheric  flow  preceded  a 
rise  in  the  barotropic  flow  in  the  ocean  by  about  3.5  days. 
This  increased  oceanic  flow  caused  an  increase  in  the  advec- 
tion  of  heat  and  oscillations  in  the  sea  surface  temperature 
about  7.6  days  later.   This  change  in  sea  surface  tempera- 
ture led  an  increase  in  the  vertical  heat  flux  between  atmos- 
phere and  ocean  and  an  adjustment  of  zonal  mean  atmospheric 
surface  temperature.   It  is  interesting  to  note  that  even 
with  the  limited  model  developed  here,  there  are  longer  per- 
iod interactions  between  the  atmosphere  and  ocean.   If  the 
model  was  modified  so  that  more  freedom  were  allowed,  then 
interactions  at  different  periods  might  occur.   There  are 
several  ways  in  which  this  model  might  be  changed  for  future 
studies.   Some  suggestions  are  listed  below. 

The  incorporation  of  continental  land  masses  and  the 
corresponding  allowances  for  a  zonal  standing  eddy  into  the 
model  is  a  necessity.   This  effect  would  modify  the  heating 
in  the  model  in  that  a  more  realistic  longitudinal  pattern 
of  vertical  heat  flux  would  emerge.   The  effect  of  this 
change  on  the  frequency  of  response  in  the  model  would  be 
interesting  to  observe.   The  incorporation  of  a  second  har- 
monic into  the  atmospheric  Fourier  formulation  would  simu- 
late the  continental  effect.   This  would  also  cause  a  shift 
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in  the  position  of  the  short  period  oscillation  or  add  a 
new  period  of  interaction.   A  study  of  an  atmospheric  model 
with  a  two-wave  structure  was  carried  out  by  Trumbower 
(1972),  showing  the  feasibility  of  including  a  second  wave 
in  the  atmospheric  formulation.   Incorporation  of  the  second 
wave  in  the  atmosphere  would  require  an  increase  in  comput- 
ing time  but  compared  to  the  total  time  required  for  the 
experiment,  it  would  be  negligible. 

A  change  in  the  sensible  heat  term  should  be  incorpora- 
ted in  the  model.   In  the  term  Hn (T     -  T.),  the  quantity 

1   sea     4  J 

H   was  held  constant.   This  quantity  should  be  divided  into 
two  parts,  one  that  is  constant  and  one  that  is  a  function 
of  the  surface  wind  speed.   This  would  cause  a  portion  of 
the  vertical  heat  flux  to  be  dependent  on  the  wind  speed 
and  allow  more  freedom  in  the  model  resulting  in  other 
forms  of  interaction. 

An  analysis  of  the  energy  budget  of  the  ocean  should  be 
incorporated  in  any  future  calculations.   The  transforma- 
tion of  energy  between  the  vertical  shear  kinetic  energy 
and  the  vertical  mean  kinetic  energy  and  the  dissipation  of 
energy  is  an  area  well  worth  studying. 

The  phase  III  experiment  should  be  carried  out  again 
with  a  sea  surface  temperature  with  an  x  dependence  inclu- 
ded.  A  zonal  variation  in  the  heating  would  be  present 
and  the  frequency  at  which  the  variables  oscillated  may  be 
changed  . 
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For  atmospheric  prediction  for  periods  shorter  than  a 
few  weeks,  an  interacting  model  may  not  be  necessary.   An 
atmosphere-only  model  with  an  accurately-observed  air-sea 
temperature  difference  may  be  sufficient.   However,  longer 
range  predictions  on  the  monthly  or  seasonal  time  scale 
should  include  some  kind  of  simultaneous  ocean  calculation. 

The  atmosphere-ocean  model  was  synchronously  integrated 
for  a  79  year  simulation  of  the  circulations  of  the  two 
media.   The  interaction  observed  at  periods  of  one  and  two 
months  provides  evidence  of  the  importance  of  the  effect  of 
the  ocean  circulation  on  the  transient  behavior  of  the  at- 
mospheric circulation,  and  vice  versa.   Other  periods  of 
interaction  surely  exist  but  the  investigation  of  all  of 
them  is  beyond  the  scope  of  this  work.   It  is  hoped  that 
the  results  presented  in  this  paper  are  sufficient  to  moti- 
vate further  research  in  the  important  realm  of  large-scale 
air-sea  interaction. 
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APPENDIX  A 
Derivations  of  Equations  (2.22)  to  (2.27) 

Equation  (2.14)  may  be  rewritten 


-^(V2iy  +lk  x  Vi|;M  •  V(VZ^M)  +  Ik  x  V^T  • 

3t 


V(V2i{;  )  +  3 


o  3x 


9  $*,                         9  T 
:  = S_r X 
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fi  p  3x     3y   4 


where 


and 


ijj   =  E(y,t)  +  A(y,t)  cos  kx  +  B(y,t)  sin  kx, 
ipT    -  F(y,t)  +  C(y,t)  cos  kx  +  D(y,t)  sin  kx , 

(Tx}4  "  P°D  Kl  U4 
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An  examination  of  quantities  in  Equations  (A4)  and  (A5) 

shows 

3iK 
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Linear  interpolation  was  used  to  obtain 

3  i     1  , 
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and  the  relationships 
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and 
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Then 
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Apply  Equations  (A2) ,  (A3)  and  (A9)  to  (A6)  to  obtain  (when 
2k  and  higher  order  terms  are  neglected) 
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(A12) 


Equation  (A12)  is  simplified  by  making  the  assumption 
that  contributions  of  the  disturbances  are  much  less  than  the 
contribution  of  the  zonal  wind  to  the  magnitude  of  the  sur- 
face wind.   The  cosine  and  sine  terms  in  (A12)  are  therefore 

neglected  and  a  new  formulation  for  IvV. I  is  obtained. 

'  '  4  ' 
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Equation  (A4)  is  now 
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When  Equations  (A2),  (A3)  and  (A16)  are  substituted  in 
Equation  (Al)  and  terms  of  wave  number  2k  and  higher  are 
neglected,  an  equation  of  the  form 

y^[G1(y,t)  +  G2(x,y,t)  cos  kx  +  G3(x,y,t)  sin  kx] 

=  F1(y,t)  +  F2(x,y,t)  cos  kx  +  F  (x,y,t)sin  kx   (A17) 

is  obtained.   Coefficients  of  the  terms  independent  of  x, 
coefficients  of  the  cosine  terms,  and  coefficients  of  the 
sine  terms  are  equated  to  yield  the  following  three 
equations  : 
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Rewrite  Equation  (2.15)  to  obtain 


3Y    (V-yz)*T   +!k    x    V^M    •    V(V2-p2)^T 

3*T  y2KQ 

2f 
o 


+  lk    x    V^T    •    V(V^M)    +    3Q    ^ 


(A21) 


0       (     3t  3t  3t  3t        \ 

ftp     T3x  3y    J2  L3x  3y    J  4 1  ' 

and    examine    the    term    u2KQ/2f     .       Substitute    for    Q    to    obtain 
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Define    a    constant,     C, ,     such    that 
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and  the  term  is  then 
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Substitute  Equation  (A3)  for  ty       and 
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is  obtained  . 
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The  terra  -g/Ap  2  [-r — ^  -  ~ ]0   is  then  examined.   Recall 
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that  can  be  broken  up  into  its  vector  components 
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The  quantity  T   may  be  written  as 
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that  leads  to 
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which    leads    to 
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Finally 
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and  can  be  written  in  its  Fourier  components  as 


A  p       2   p  2  gy2 


3'C 

3y2 


)cos  kx-(Dk  -  )sin  kx] . 

3y2 


(A35) 
Repeat  the  procedure  previously  used  on  Equation  (Al) 
on  Equation  (A21)  to  obtain  three  equations.   These  are 


3  r32C    .  2„    2ot   ,  r3E  3D  .  3F  32B    3E  _  „  7 
^[—  -  k  C  -U2C]=    k[-y  —  +__-_  D(k^  +  U2) 


3y 

3F   ..  2    2S 
-  ^—  B(k  -  y  ) 
dy 


33E 


3y 


3y 

3  3F 


D  - 


3y 


B]  -  3   Dk 

3  O 


p2k 


2f 


A  P  Hi(-  -TTC>  +  ~fc   ^"M^MM  -   ck2> 

o  r  3y 


+ 


■4 


|\V.  |k2(2C-A)-  — [  |VV.  |  (2-|£  -  |^ 
Ap    D  | '  4  '  3y  '  4  '   3y    3y 


».] 


3   r3  D    ,  2  ~    2ni  i  r  3  E  _  / ,  2  ,   2\,   3  F  .  , .  2       2. 

j^[ k':D-yzD]=k[^-C(k^  +  y^)+  —  A(k   -  u  ) 


3y 


3E32C    3F32A,33F 


dy    3y2  8y    3y2    +    3y3 


A    +    C]    +    3       Ck 

r,  3  O 

3y 


o  r  r  gy* 

+  "aV  pcd  [f4l  k2(2^-B)-  |yfKV4|(2|f  -  |f)j] 


(A36) 


(A37) 


2  0.1 


and 


at 


By 


/  J  dy 


3  r  B  2  D 


A  -  B- 


82C 


By 


+  c 


B2B 

3y2 


-  D 


a'A,    .   .2  9 


By 


J  +  y 

2        By 


By' 
[BC  -  AD] 


y2Kg 

2f  Ap  "l^sea 
o 


2f 


H.  [T 


^(F  +  Cl>] 


3y 


A  p 


pcj-rlvjc3 


DS 


4  '  By    By 


2^)]. 


(A38) 


This  completes  the  derivation  of  the  equations  for  pre- 
dicting new  values  of  the  Fourier  coefficients. 
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APPENDIX  B 


Derivation  of  w  and  v..  Equations, 


Applied  at  level  two,  Equation  (2.12)  is  written  as 


3    ^  "  ^1 


>  +  f-»2"  "  f^?  V  <B1> 
o  o 


A  p   y  "  '"  A    "y2  '  "   A  p 

Substitute    from    Equation     (2.13)    which    yields 

g  A   p  K 


d 

3t    rT 


i+lkx    V*M    -    V^         ff-    •   w2    =    2f~   Q2     ' 


M 


2f 


(B2) 


Solve    for    co„    which    give. 


2f 


3^T       3^M  3ipT       a*M  3tpT 


W2         aip     \3t  3y       3x       '     3x       3y 


IT*2 


(B3) 


Substitute  Equation  (2.16)  and  (2.17)  into  (B3)  and 
neglect  terms  of  wave  number  2k  and  higher  which  gives 

2fo   (lM  ,_3C    ^3AN    .  ,  3B    A3D,, 

w0  =  — -r — {   -r-[k(B-r—  -  D^—)  +  k(C-r A-5— )  ] 

2    a  A  p  l  2     3y     3y         3y     3y 

,  3F    KQ2   r3C    .  ,9E_    3Fn.  .      . 
+  9T  ~  7f~  +[3l  "  k(37D  ~  37B)]  COS  kx 

o  J  J 

+  [|t-  +  k(~C  -  |^A)]  sin  kx 
d  t       3y     3y 

or 


03„  =  to  +  (00  )  cos  kx  +  C^O  sin  kx . 


(B4) 


The  quantity  oo„  is  composed  of  a  mean  component  and  a 
disturbance  component.   The  mean  component  may  be  written 

77   -        £o       M  ^p9c         n9^    j.    1   /^9B  A3D, 

w         nAn     lk^B^777   "    D"5~)    +    MCtt—   -    A^— ) 
0  A  p  dy  dy  3y  3y 

t 

3F         K<^2 
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and  the  disturbance  part  written  as  two  components 
2£o   r3C    W3E_    3FRx1 


wc  = 


and 


w 


J) 


2f 

o  r  3D  .  /3Er,  3F  ..  , 

a IT-  +  k(-r—  C  -  -r— A)  ]  . 

a  A  p  l  3t      3y  3y 


(B6) 


(B7) 


The  mean  meridional  motion  component  can  be  derived 
from  the  zonally  averaged  continuity  equation.  This  is 
writ  ten 

(B8) 


V  •  W     +  ~    =  0. 
3p 


Applied  at  level  one  this  yields 


3v. 
37 


0), 


A  p 


(B9) 


and  this  is  solved  numerically  to  obtain  v       at  any  point  in 
the  model,  after  obtaining  go-  from  (B5). 
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APPENDIX  C 
Derivation  of  Energy  and  Energy  Transformation  Equations. 

To  derive  the  energy  relations hips,  the  A,  B,  C,  D,  E 
and  F  equations  [(2.22)  to  (2.27)]  were  multiplied  by  A 
through  F,  respectively.   After  further  mathematical  mani- 
pulation, the  energy  equations  were  obtained. 

Multiply  the  A  equation  by  A  to  obtain 


A3  r32A 
3y 


.,2,    .  tr3E  32B  x  3F  32D    33E,, 
-  Ak    -  1cA[-a—  +  -r—  -         — B 

3y  ^2    3y 


3y 


3y2    3y3 


33F 


3y 


D  -  k2(~B  -  |^D)J  -  3  ABk 
3         dy      dy         o 


PC,     |\V4|  k2(A-2C)  +|-IKV4|(2|£.  M)] 


(CI) 


and  the  B  equation  by  B,  and  add  to  Equation  (CI)  to  get 


3  r32A 


kJtl 


3y 


-    Ak  ]+  Bt— 

2  d  t 


3  r32B 


3y; 


,2,      A3E  32B 
ik  ]  =  k[A-r—  

3y  3y2 


„3E  d*k  .  3F  32D 

dy  „  2      3y  ^  2 
3y       J     3y 


3F  32C    AT,33F 
3y  3y2      3y3 


+  BC 


33F 


3y 


-  k2(|^  DA  -  ~    CB) ] 

3      3y      3y 


A?  PCD 


|\V,  |  k2[(A2-2AC)  +  (B2  -2BC)] 


3C    3A 


3D    3B- 


+  A^l  K  I  (2^  "  ie>  J  +  »f-[  [W,.  |  (2^  -  |£)  i  . 


3y  '  4  '   3y    3y 


3y 


9y    3y 


(C2) 
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Treat  the  C  and  D  equations  similarly  and  add  to  obtain 


3tla  2 

3y 


-  Ck2  -  Cy2]  +  d|-[^ 

dt    By2 


-  Dk   -  Dy2]  - 


kI|l(Ci^  -  D^)  +  ^(C^  -  D^) 

oV    ,,   2        ^2        dY^2        •>   2 

J        dy      3y       '   dy      3y 


33F 


3y 


(BC  -  AD)  -  (k2-  y2)  l^(BC  -  AD) 

3  dy 


+  f?  PCD  JKI  k2[(2C2-AC)  +  (2D2-BB)]-c|7[|W4|(2|£-  |A,  ] 


•I>|-[|\V.  |(2f£ 


3y;j  (       a  P 


4(pic),(£J)9[C(^ 

2      p    2         9y2 


-Ck2) 


32D 


+    D(— -    -    Dk    )  ]  . 

3y2 


(C3) 


Combine    Equations     (C2)     and     (C3)     and    integrate    over    the 


region    to    obtain 

y    , 


/ 


AtttI -Ak    ]+B^—[ -Bk    ]+(hr-[ -Ck    -Cy     ]  + 

dy  dy  3y 


D|_^   -Dk2-Dy2j       dy= 

9y2  ) 


y>n 


(A 


32B       1,32A 

—  is 

3y2  3y2 


+    C 


32D 


dy 


D^) 

3y2 

+   |Z(AilD    „    D 

^       3y2 


23F 


+(B-D) 


3y2         3y2         3y2  dy 

4 


[(A-C)  flw,.  |k2(A-2C)+       rllW,  I  C2-| 


3y 


3y         3y 


^1^(3^)^11^1(2^-^)])] 
+    -ih    4<PtOo(%oIC(—    -Ck2)+D(^    -Dk2)]  1   dy. 


2  v    p'  2 


3y 


ay 


I 


(C4) 
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Perform  the  multiplications  and  apply  the  boundary  con- 


ditions to  acquire 

T 

d 


T-/ 


,9A,  2^,dS>.  2J_,DCk2,  ,c)D,.2,  ,  2/A2,_2x.2,_2.  ,  2,-2.  _2X. 
^3"" )  +(T~ )  +(T~)  +(^t)  +  k  CA  +B  +  C  +D  )+u  (C  +D  )  dy 

dy     dy     dy     dy 


/ 


|[.  |E(A8lB  .  BifA  +  cliD  _  ^    |F(AiiB  _  DilA 

2  DY       r>    2  n.    2  rs    2  r^2  <?  >'       *N    2  ^2 

dy      dy      dy      dy         dy      dy 


-  B-^  +  c^-^-)-y2|^(BC-AD)]-  -^  pC 


^2       ^2       dy 

dy      8y       J 


Ap 


D   [<A"C>( 


IV 4  |k2(A-2C) 


+  f7[KI(2lf  "  I7)]  +(»-«)( Klk2(B-2D)+  f-[|wj(2M 


3B 


>)] 


32C   „,  2X  ,  ^,32D   _  2 


v)])|    T^  A(pic)  (%  [C(-    -Ck")+D(-    -Dk")]>dy 
Ap        2   P  2    9y2  9y2  j 


(C5) 
Apply  similar  techniques  to  the  E  and  F  equations  to 
ob  tain 


t/ 


fe/  i<S>'.#!-»''!i« 


k[|E(AiiB    .    B3lA    +    c3iD 
2    ^       3y2  3y2  3y2 


„    Di!£)+    3F      £D    _t3lA    _B3!c    +cl!l         2|F  D)] 

3y2         ^      3y2         3y2         3y2         3y2  3y 

r    o 


^-^-^w.^,^11" 


(C6) 
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Define 
P, 

P 

< 

K 

K 


D  . 
ei 


es 


zi 


zonal  potential  energy, 
eddy  potential  energy, 
zonal  kinetic  energy, 
eddy  kinetic  energy, 


conversion  of  P   to  P  , 

z      e 

conversion  of  P   to  K  , 

e      e 

conversion  of  K   to  K  , 

e      z 

conversion  of  K   to  P  , 

z      z 

generation  of  P   by  diabatic  processes, 

dissipation  of  K   by  internal  friction, 

r  e  J 

dissipation  of  K   by  surface  friction, 

dissipation  of  K   by  internal  friction 
r  z 


and 


zs 


dissipation  of  K   by  surface  friction. 

z 


From  Phillips  (1956)  and  the  basic  definitions  for  i|j 
and  ty.,    used  in  this  work  (Equation  2.13)  equations  for  the 
potential  and  kinetic  energy  components  can  be  written  as 


and 


P   dy  - 


Ady 


-h 


fh2 


CC2  +  D2)  dy, 


1  /\,dE,.2,  ,ol\2.. 

2  /  [  (Tv}     ^ 


1  /"   dA  2   tfB 


)2+(|c   c|£)2 

dy    oy 


2   2   2   2   2 
+k  (A  +B  +C  +D  ) ]dy. 


(C7) 

(C8) 
(C9) 


(CIO) 
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The  right  hand  sides  of  Equations  (C5)  and  (C6)  yield  ex- 
pressions for  energy  generation,  conversion  and  dissipation 
as  follows : 


D 


D 


es 


D 


zi 


2  Apf    1   sea 


2f 


^  (F  +  CX)]F  \   dy, 


[\    y2(BC-  AD)       dy, 


II  ■ 
3y 


Ceil) 

(C12) 


(kr3E  ,A32B    _3 
MtT  (a "  -  B- 


A 


\ 


2  L3y 


9y 


3F  ,.32D 
+  ^—  (A 

8y    3y2 


9y^ 

32A 

ay2 


+  c 


32D 


^      2 

3y 
32C 

9y2 


9y2 

+  C^)]  (dy,  (C13) 


3y 


Ap"  ^P^)2(£f)2[C^  -Ck2)+D(-^ 


-Dlt)Udy,(C14) 


Sy 


3y 


£_ 


PC 


Ap  K~D 


T(A-C)  (|vV4|k2(A-2C)+  |^]| 


3yJ  n"4 
3D 


3y    By'^' 


+  (B-C)(^Kv4|k2(B-2D)+  |_[|vv4|(2|^  -  f|)])]dy,(C15) 


,32F 


f-  A(PK)  (£|)   [F-  -]  I   dy, 
Ap       2   P  2    3y2   ( 


(C16) 


and 

D 
z  s 


■/ 


fc  ps 


(E-F)  f-lKvJ(2f 


y 


9E 

3y 


)] 


dy.   (C17) 


The  quas i-geo s tr oph ic  first  law  of  thermodynamics, 
Equation  (2.8),  can  be  manipulated  in  a  similar  manner  to 


derive  expressions  for  C   and  C  .   These  are 

e       z 


and 
C 


2f( 

a7 


CO  F  dy  , 


I  A?  (WCC  +  %D)]  d?' 


(C18) 


(C19) 
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APPENDIX  D 

The  difference  form  of  the  equation  for  the  prediction 
of  the  vertical  shear  current,  u',  analogous  to  Equation 
(3.8)  is 


,  (£+1)      , 
u '        -  u  .  . , 
ijk ijk 

2At 


(£-1) 


,£     1   ^   r*    , 
'ijk    H  2Li  ijk  ^ 


Az) 


K=l 


+  AV,  f  .   (av!  .,  (£+1)  +  (l-a)v!  .,  (A  1}) 
k  j      ijk  ijk 


(Dl) 


where 


ijk    (Az) 


AA 
(I-)..,   -  -. —  [w..+w     +w.    +w      .   )   1 
1  xjk    4    ij   i+lj   ij+1   i+l,j-l  k-- 


(u).        1     -      (W.    .+W.,_    .+W.    .       ,+W...        -a-l)^1  (U).,l].. 

k--  ij       i+lj       ij+1       l+l, j+l'k+j  "k+j'ij 


7-  ^(Az>k  7[<pi+i3+*i+i.j+i>-<pij+piJ+i>]k 


+    ,A    x       )  (Az).        Au[(u. ,n  .     -    u.  . )-(u.  .    -    u.     ,  .) 
(Az)       )  yk         M         i+lj  ij  ij  i-lj 


+    (u.  -    u..)-(u..     -    u..        )] 

ij+1  ij  ij  ij-1       k 


k(uk-l~Uk)ii 

+  aa      r J    - 

l        Uz),      1 
K    2 


(Az)k+|  ] 


(D2) 


and 


AV.     =    Ax    Ay(Az) 
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and 


AA  -  Ax  Ay 


and 


(V1Jk 


^(u  .  ,+u  .  ,  -  .  ) , 
2   ij   l+l j  k 


Ay(Az) 
8 


(u  .  ,  -  .+  u  .  n  .  ,  ) 

i+lj    i+lj-1 


+  (u.  .  +  u.  .  .)  +  (u.,.,,.  +  u.   .) 
ij     ij-l       i+lj+1     i+lj 


+  (u.     +  u. .) J  ] 
±2+1  ij  /k 


-  [f^i-u  +  -ij) 


Ay(Az) 
8 


(u.  .  +  u.  .  .) 
ij     ij-l 


+  (u.  ,  .+  u.  ,  .  ..  )  +  (u  .  .  ,  -  +  u..)  +  (u.  -  .  - 
i-lj    i-lj-1'    ij+1     xj     x-lj+1 


+"i-lj>  k3 


+  [T(u. .+  u. . ,n) 
2   ij    xj+1 


Ax(Az) 


(vij+l  +  ^-Ij+l5 


+  (v..+  v.  , .)+(v. ,- .+  v. .)+(v. .- .  .+  v,  .  . )   .  ] 

ij    i-lj     i+lj    ij     i+lj-1    xj-1  /k 


+  -(u. .  .+  u. .) 
2   ij-1    ij 


Ax(Az) 
_ 


(v  .  .  +  V  .   ,  .  ) 

IJ    l-lj 


+(v..  -+v.  -.  ,)+(v. .. ,+v. .)+(v. .-.  ,+v,.   )    ]. 

xj-1   x-lj-1     x+lj   ij     x+lj-1   ij-1  /  k 

(D3) 
The  difference  form  of  the  equation  of  the  meridional 
component  of  the  vertical  shear  current,  v',  analogous  to 
Equation  (3.9)  is 


v,  (A+l)   ,  (£-1) 

ijk ijk 

2At 


=  R 


I  1   fe  _£ 

ijk    H  2L,        ij 


jk  (Az)k 


k  =  l 
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■AV.  f  .[a(u')  ^1}    +    (1+a)    u'    ^"1)] 
k    j  ljk  iJK 


(DA) 


where 
xjk 


1         (  AA 

-7-7 — r —    J  -  (I,. )  .  .  ,       -    -: — [  (w  .  ,+W.  ,  -  ,+W.  .  ,  ,+W.  ,  -  .  ,  n  ).       1 
(Az)         )  I     Xjk  4  ij         X  +  lj         Xj+1        X+lj+1     k-y 


(v)k_l    -(-i:j+-i+lj+wij  +  14-wi  +  1.+1)k+l     (v)k+l]±j 


-    —   Ax(Az).     y[  (P.  ..,+P.^ ux1  )    -     (P,  .+P4X1  .)],   j 
P„  k    2lv    13+I       i+lj+1'  xj       x+lj/Jk/ 


(Az) 


(Az).       A_,[(v. ,n .-v. .)    -     (v..    -    v.     .,) 
k         M         x  +  lj        xj  xj  1_lj 


+     (v.,  ..-    v.  .)     -     (v.  .     -    V.  .     .  )  ]. 
ij+1         xj'  xj  xj-l/Jk 


+    AA[ 


k(vk-i-  Vn     k(vk  -  vk+i>ii 


(4«)k  1 


and 


(IJ  .  .,      =     [~(v.  ,+v. .-  .), 
2    ljk       )    2       xj        l+l J    k 


(Az)k4 


Ay(Az) 


■] 


(u  .  ,  1  .     +    u  .         .     .  ) 
x  +  lj  x+lj-1 


(D5) 


+     (u..+    u..        )     +     (u.  ,-.,.)     +     (u, . ,,+    u..)        ] 
xj  xj-1  x+lj+1  xj  +  1         xj       k 


[T(v.      .   .+    V. .) 
L2       x-lj  xj 


Ay(Az) 


(u . .+    u . .     . )+(u .     - .+u .     . .     , ) 
xj  xj-1  x-lj       x-lj-1 


+     (u,  .  .-+    u.   .)     +     (u.      -  .  ,-  +    u.      ..)).] 

xj+1  xj  x-lj  +  1  x-lj'/k"1 


+  i^u^ij+i) 


Ax(Az) 


(vij+l+  vi-lJ+l> 


+     (v. .+    v.     _  .)  +  (v.  ,  .  .+    v.  .  )  +  (v.  ,  .  .     -  +    V.  .     .  )     .   ] 
xj  x-lj  x+lj  xj  x+lj-1         xj-1     /k 


+  [o"(v.   .      -.      +v.    .) 

2       xj-1  xj 


Ax(Az) 


k 


(  (v.  .+    v.     . .) 
\       *0  x-lj' 


+  (v..     ,+v.     ,.     ,)  +  (v.  ,  n  ,+v.  .)  +  (v._Ll  .     ,+v..     _)).]}.  (D6) 
xj-1       i-lj-1'  i+lj       ijy  x+lj-1       lj-l'/k"1 
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At  each  time  level  £ ,    Equations  (Dl)  and  (D4)  are  solved 


£+1 
together  for  the  unknown  shear  currents,  u ! ..       and 

ij  k 

,  £+1 
v' 
ijk 

The  difference  form  of  the  vorticity  equation  analogous 
to  Equation  (3.12)  is 


£+1 


l-l     rl 


/ 


ij         ij       2\   3  ij    3  ij-1 


I   (<Vi-lj  +(Vi-l,:-l)  -  r(<Vlj  +  <Vi-lj 

+  I  /  (I  )    +(I  )      \  _  3  (  1+1J 1"1J)] 

2  yvx4;lj-l  ^,L4;i-lj-li    PoV    2Ax       ;J 


where 


km 


<13>ij  -IE    i 

k  =  l 


(D7) 


-   S-(I1)£.,  +  (Az),  A M[(u. ,_ .-u. .) 
y   J     1  ijk       k   M    i+lj   ij 


(u..-  u.  1,)  +  (u..,1-  u..)-(u..-  u..  _)]<*  1} 


+  AA  ^Vj-vn .  t(vyii]<wi) 


k    2 


(Az)k+i 

k    2 


(D8) 


(I,) 


,(-I0)7.1.+  (Az).     AMI(v.,1  .    -    v.  .) 

k       M  l+l j  ij 


'ij         H      2     Ax  r    J"2/iJk 
k=l 


-     (v..-v.     .  .)  +  (v.  ,.,-v.  .)-(v.  .-v.  .     ,)]$Z    1) 
ij       i-lj  ij+1       ij  ij       ij-l/Jk 


k(v,      ,  -    v.  )  .  .  k(v.  -    v.        )  .  .     0     .. 

aat         k-1         k    ij  k         k+1    ij,£-l 

k 


(Az)       1 
k    2 


(Az)k+| 


and 


(D9) 


ij       2iri+lj+l    ri-lj+l    ri+lj-l    ^i-lj-1  ij 


(D10) 


obtain  \p 


Equation  (D7)  is  solved  by  relaxation  techniques  to 

£+1 


ij 
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The  difference  form  of  the  thermodynamic  equation  analc 
go us  to  Equation  (3.16)  is 


<T..  <1+1>-  T..  C*"1)) 

2At 


A\   =    -    "Vijk   +(wij+  Wi+lj 


+    Wij+1    +    Wi+lj+l>k-f    (T)k4   -     ^ij+'i+lj^ij+l 
+    -i+lj+l)k+|    Wk+|J*    +     <Az>k    AHf(Ti  +  lj     "    T±J> 


-     (T..     -    T.     ,.)+(T..,.-    T..)-(T..     -    T..        )] 
13  i-lj  ij+l  ij  13  iJ-1       k 


(£-1) 


k(Tk-l~    V  k(Tk    "    Tk-1}  (£-1) 

(Az)n    1  (Az)nj_l  Jij 


k- 


k+2 


(Az) 


-    1 


+ 
where 


(I,).,,      ^v[(T..+    T.,..).      •     (u..+    u..     .). 
5    ljk       /4         ij  i+Ij    k  ij  ij-1'k 


(DID 


-     (T.        .+    T. .).  •     (u.     . .+    u.     , .     n ). 
1-I3          13 'k  i-lj  i-lj-l'k 

+    (T.  .+    T.  )  •     (v.  .+    v.     ,  .). 

ij  ij+l    k  ij  1-I3    k 


and 
(Az) 


At  v~~' C 

vious ly . 


-     (T.  .        +    T.     )        •     (v.  .     .+    v.     .  .     .  ).   ] 
ij-1  13     k  13-I  i-l3-l/kJ 


(6T)        is    the    convective    adjustment    discussed    pre- 
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A  two-level,  quasi-geos trophic  atmospheric  model  on 
latitude  B-Plane  with  a  single  wave  in  the  x-direction  was  coupled 
with  a  three  level  primitive  equation  ocean  model  through  the 
exchange  of  sensible  heat  and  momentum  to  investigate  the  eftect 
the  ocean  circulation  on  transients  in  the  atmospheric  circula- 
and  vice  versa.   A  long-term  (79  years)  simulation  of  the 
atmosphere-ocean  system  was  carried  out.   The  experiment  was  con- 
ducted in  three  phases  that  consisted  of  integrating  the  combined, 
the  ocean-only  and  the  atmosphere-only  model  equations. 

The  combined  model  produced  realistic  patterns  of  both  atmos 
pheric  and  oceanic  circulations.   Energy  studies  of  the  atmosphere 
were  completed  with  the  results  comparing  favorably  with  the 
existing  knowledge  of  the  energy  balance  in  the  real  atmosphere. 

Data  were  retained  at  selected  points  on  a  daily  basis  and 
spectrally  analyzed  for  a  simulated  six  year  period.   Transients 
with  periods  on  the  order  of  one  and  two  months  appeared  in  the 
spectra  of  both  atmospheric  and  oceanic  variables  from  the  combi 
model.   These  transients  were  not  observed  in  the  ocean-only 
atmosphere-only  phases  of  the  experiment,  and  are 
interaction  occurring  in  the  combined  model. 
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